/****************************************************
 * File: pick.c
 *
 * Author: Fred Engler
 *
 * Description: Implementation of Dijkstra's shortest
 * paths alogorithm for a sparse matrix, and
 * collection/printing of stats for the current run.
 * Originally part of standalone PickMTP package.
 * Modified to read from 'pairs' struct instead of from
 * file, and to assign different weights based of the pair
 * type (fingerprint or BSS-draft).
 ****************************************************/

/* cari 3 feb 05 - change all BANDLEN to Proj.avgbandsize */

#include "pick.h"
#include "fpp.h"
#include "glib_port.h"
#include "clam.h"
#include "mtp.h"
#include <string.h>
#include <stdio.h>

extern void mtpParameters();
extern int *pick_cnt; /* cari March 06 - need for write_stats */
extern struct picklist *picks;
extern int *pair_cnt;
extern struct pairlist *pairs;
extern int maxShared;
extern int prefer_large;
extern int run_on_single_ctg;
extern int selected_ctg;

extern struct pair* find_pair(int ctg, char *c1, char *c2);
extern void add_to_picks();
extern int loadCz(struct tmpCz *Cx, int cin);
extern int GetSize();
extern void init_ids();

char* clone_idx2name(int i);

typedef struct
{
    guint left;
    guint right;
    EdgeWeight score;
    guint subgraph_idx;
} ClonePath_t;

typedef struct
{
    guint sm_index;
    gboolean in_subgraph;
    gboolean in_this_subgraph;
    gint left;
} gap_clone_t;

typedef struct
{
    struct subgraph *sgf;
    gint *parent;
    guint sm_left;
    guint sm_right;
    guint sp_left;
    guint sp_right;
    gboolean left_pref;
    gboolean right_pref;
    EdgeWeight score;
} gap_path_t;

BOOL fppFind(Array a, char *s, int *ip, BOOL (* order)());
BOOL cloneOrder(void *s, void *b);
static void add_to_results(int ctgID, struct subgraph *sgf,
                           int *parent, guint left, guint right);
static void free_mandatory_subgraph(struct subgraph *sbgrph,
                                    gpointer reserved);
static gint compare_left(const CLONE **c1, const CLONE **c2);
static void add_mandatory_to_results(guint ctgID, struct subgraph *sgf);
static void span_gap(guint ctgID, GArray *left_ends, GArray *right_ends);
static GArray* find_gap_clones(guint ctgID, GArray *left_ends,
                               GArray *right_ends);
static gint compare_gap_clones(const gap_clone_t *gc1, const gap_clone_t *gc2);
static GArray* find_gap_subgraphs(guint ctgID, GArray *gap_clones,
                                  GArray *left_ends);
static void scan_gap_subgraph(const struct vertex *sm, GArray *gap_clones,
                              gap_clone_t *start, GArray *subgraphs);
static GArray* find_spanning_paths(guint ctgID, GArray *gap_subgraphs,
                                   GArray *left_ends, GArray *right_ends);
static gint paths_compare(const guint *idx1, const guint *idx2,
                          GArray *gap_paths);
static guint best_path(guint ctgID, GSList *candidates, GArray *gap_paths,
                       gboolean prefer_left);
static guint leftmost_path(guint ctgID, GSList *candidates, GArray *gap_paths);
static gboolean is_past_ends(struct vertex *sm, gap_path_t *path,
                                  GArray *right_ends);
static void trim_and_add_paths(guint ctgID, GArray *spanning_paths);
static gboolean path_ok(guint ctgID, gap_path_t *gp);
static gint compare_subgraphs_left(const struct subgraph *s1,
                                   const struct subgraph *s2,
                                   struct vertex *sm);
static void join_expressways(guint n_ctgs);
static gint compare_expressways(const struct expdata *e1,
                                const struct expdata *e2);
static void shortest_path(guint ctgID, EdgeWeight *dist, GArray *left_ends,
                          GArray *right_ends, gap_path_t *gp);

void pick_mtp(char *outfile, int write_xl_file);
void pick_exp(int output_fpc);
void load_sparse_matrix(void);
extern int use_mandatory(CLONE* clone);

int max_weight;
int min_size_pen;
int max_size_pen;
int min_olap_pen;
int max_olap_pen;

struct mtpdata *mtp[MAXCTGS];
/* CAS =  uses too much memory
char clonekey[MAXCTGS][MAXCTGSCLONE][16]
*/
int MaxCtgClones;
int *clonekey[MAXCTGS];
int num_clones[MAXCTGS];

int ctgkey[MAXCTGS][3];
int num_ctgs=0;

struct vertex *SMctg[MAXCTGS];

int exp_num=0;

int total_map_size=0;

extern gint shared_weight, unacc_weight, spanmis_weight;

static GArray *clones_ctg = NULL;

#define UDIV_RND(n, d) ((2 * (n) + (d)) / (2 * (d)))


/**************************************
        pick_exp_and_mtp
CAS 19oct03
***************************************/
void pick_exp_and_mtp()
{
    int i;
    GArray *null_array = NULL;

   max_weight = maxShared * Proj.avgbandsize;
   min_size_pen=INT_MAX;
   max_size_pen=INT_MIN;
   min_olap_pen=INT_MAX;
   max_olap_pen=INT_MIN;

/*  for (i=1; i<=max_contig; i++) {
    if (MaxCtgClones < contigs[i].count)
        MaxCtgClones = contigs[i].count;

    clonekey[i]  = malloc(contigs[i].count * sizeof(int)); // WMN 10/10/06 
  }   WMN This block should be taken care of in init_ids */

  if (MaxCtgClones > MAXCTGCLONES) {
     printf("Contig has maximum of %d clones\n",MaxCtgClones);
     printf("Change MAXCTGCLONES in mtp/mtp.h\n");
     printf("or email cari@genome.arizona.edu\n");
     messcrash("Too many clones per contig");
  }
  init_ids();
  clones_ctg = g_array_new(FALSE, FALSE, sizeof(GArray*));

  for (i = 0; i <= max_contig; ++i) {
    GArray *clone_array;
    struct contig *ctg;

    g_array_insert_val(clones_ctg, i, null_array);
    if ( i == 0 || contigs[i].count == 0)
        continue;
    clone_array = g_array_new(FALSE, FALSE, sizeof(CLONE*));
    g_array_insert_val(clones_ctg, i, clone_array);
    find_Clone(i);
    orderclones(i);
    for (ctg = root; ctg != NULL; ctg = ctg->new) {
        CLONE *cln = arrp(acedata, ctg->next, CLONE);
        g_array_append_val(clone_array, cln);
    }
  }



  printf("\n************ Starting PickMTP ************\n");
  pick_exp(0);
  pick_mtp("stdout",0);
  printf("\n");
  mtpParameters(stdout);
  printf("************ Finished PickMTP ************\n");

  /*for (i=0; i<=max_contig; i++) {  WMN */
  for (i=0; i<=num_ctgs; i++) {
      if (clonekey[i] != NULL) {
          free(clonekey[i]);
          clonekey[i] = 0;
      }   
  }
  if (clones_ctg != NULL) {
      guint j;
      for (j = 0; j < clones_ctg->len; ++j) {
          GArray *ary = g_array_index(clones_ctg, GArray*, j);
          if (ary != NULL) {
              g_array_free(ary, TRUE);
              ary = 0;
          }
      }
      g_array_free(clones_ctg, TRUE);
      clones_ctg = NULL;
  }
}
/********************************
 Input is a ctg and clone name.  Output is either its existing id or a newly created one
*********************************/

int get_clone_id(int ctgID, char *c){
   int i, index;

   if(!fppFind(acedata,c,&index, cloneOrder)) {
       printf("FPC ERROR could not find clone %s \n",c);
       exit(0);
   }

   for(i=0;i<num_clones[ctgID];i++){
      if(index == clonekey[ctgID][i])
         return i;
   }
   i = num_clones[ctgID];
   clonekey[ctgID][i] = index;
   num_clones[ctgID]++;
   if ( num_clones[ctgID] > contigs[ctgkey[ctgID][0]].count)
    {
        g_assert(0);
    } /* WMN 10/10/06 */

   return (num_clones[ctgID]-1);
}

/********************************
 Input is a ctg number.  Output is either its existing id or a newly created one
*********************************/
int get_ctg_id(int ctg){
   int i;

   if (num_ctgs == MAXCTGS) {
      printf("Too many contigs for MTP (max %d)\n",MAXCTGS);fflush(0);
      exit(0);
   }
   for(i=0;i<num_ctgs;i++){
      if(ctg==ctgkey[i][0])
         return i;
   }
   ctgkey[num_ctgs][0]=ctg;
   ctgkey[num_ctgs][1]=contigs[ctg].left;
   ctgkey[num_ctgs][2]=contigs[ctg].right;
   clonekey[num_ctgs] = malloc(contigs[ctg].count*sizeof(int));
   if (clonekey[num_ctgs] == NULL) messcrash("Cannot alloc memory for clonekey.");
   num_clones[num_ctgs] = 0;
   num_ctgs++;
   return (num_ctgs-1);
}

/********************************
 Insert edge from vertex v1 to vertex v2 into sparse matrix for ctgID.
*********************************/
void
insert_sparse(int ctgID, int v1, int l1, int r1, int v2, int l2, int r2,
              long bssOlap, int overlap_score, enum pairTypes pair_type,
              int sizeL, int sizeR)
{
    struct edge *p;

    if (v1 >= MaxCtgClones || v2 >= MaxCtgClones) {
        fprintf(stderr, "\n** Too many clones in contig. Change in "
                "pick.h and recompile.\n");
        fprintf(stderr, "or email cari@genome.arizona.edu\n");
        exit(0);
    }
    if(SMctg[ctgID][v1].count == -2) SMctg[ctgID][v1].count = 0;
    if(SMctg[ctgID][v2].count == -2) SMctg[ctgID][v2].count = 0;
    p = (struct edge *)malloc(sizeof(struct edge));
    if (p==NULL) messcrash("Out of memory on insert_sparse\n");

    p->num=v2;
    p->bssOlap=bssOlap;
    p->overlap_score = overlap_score;
    p->fpcOlap=r1-l2+1;
    p->pair_type=pair_type;
    p->sizeL=sizeL;
    p->sizeR=sizeR;
    p->next=SMctg[ctgID][v1].link;
    SMctg[ctgID][v1].link=p;
    SMctg[ctgID][v2].count++;

    SMctg[ctgID][v1].l=l1;
    SMctg[ctgID][v1].r=r1;
    SMctg[ctgID][v2].l=l2;
    SMctg[ctgID][v2].r=r2;
}

void
insert_sparse_vertex(int ctgID, int vtx, int x, int y)
{
    if (vtx >= MaxCtgClones) {
        fprintf(stderr,
                "\n** Too many clones in contig. "
                "Change in pick.h and recompile.\n");
        fprintf(stderr, "or email cari@genome.arizona.edu\n");
        exit(0);
    }
    if (SMctg[ctgID][vtx].count == -2)
        SMctg[ctgID][vtx].count = 0;

    SMctg[ctgID][vtx].l = x;
    SMctg[ctgID][vtx].r = y;
}

/********************************
Load sparse matrix from the data in the pairs structure.
Get clone, ctg ids, and insert edges into sparse matrix.
*********************************/
void
load_sparse_matrix(void)
{
    int ctgID, c1ID, c2ID;
    int i, j, k;

    total_map_size = 0;
    if (pair_cnt != NULL) {
        for (ctgID = 0; ctgID < num_ctgs; ctgID++) {
            i = ctgkey[ctgID][0];
        /*for (i = 1; i <= max_contig; i++) {   WMN */
            GArray *clone_array;
            if (run_on_single_ctg && i != selected_ctg)
                continue;
            /*ctgID = get_ctg_id(i);*/
            assert(ctgID ==  get_ctg_id(i)); /*WMN*/
            SMctg[ctgID] = g_new(struct vertex, MaxCtgClones);

            for (k = 0; k < MaxCtgClones; k++) {
                SMctg[ctgID][k].count = -2;
                SMctg[ctgID][k].link = NULL;
            }

            clone_array = g_array_index(clones_ctg, GArray*, i);
            if (clone_array != NULL) {
                for (k = 0; k < clone_array->len; ++k) {
                    CLONE *cln = g_array_index(clone_array, CLONE*, k);
                    int cln_id = get_clone_id(ctgID, cln->clone);
                    insert_sparse_vertex(ctgID, cln_id, cln->x, cln->y);
                }

                total_map_size += contigs[i].right - contigs[i].left + 1;
                if (pair_cnt[i] > 0) {
                    printf("\rBuilding graph for ctg%d...",i);
                    fflush(stdout);
                    for (j = 0; j < pair_cnt[i]; j++) {
                        c1ID = get_clone_id(ctgID, pairs[i].p[j].c1);
                        c2ID = get_clone_id(ctgID, pairs[i].p[j].c2);
#ifdef MTP_DEBUG
                        //fprintf(stderr,"SPARSE PAIR:%d %d\n",c1ID,c2ID);
#endif
                        insert_sparse(ctgID, c1ID,
                                      pairs[i].p[j].x1, pairs[i].p[j].y1,
                                      c2ID, pairs[i].p[j].x2, pairs[i].p[j].y2,
                                      pairs[i].p[j].olap,
                                      pairs[i].p[j].weight_score, pairs[i].p[j].pair_type,
                                      pairs[i].p[j].c1size,
                                      pairs[i].p[j].c2size);
                    }
                }
            }
        }
    }
    printf("\rBuilding graphs completed.         \n");
}

/********************************
 For a given ctg and first vertex, find the final vertex, total number of vertices,
 and all vertices in the subgraph reachable from first.
*********************************/
void scan_subgraph(int ctgID, int first, struct subgraph *sub_info){
   struct vertex *local_a;      /*A local copy of the sparse matrix for ctgID
                                  Only the vertices are copied -- the edges are not altered;
                                  hence we use a pointer to them.*/
   int start, i, j, k;
   struct edge *p;
   int last;             /*Vertex with most extreme right end value*/
   int max_right=0;      /*Highest right end value*/

   local_a = (struct vertex *)malloc(MaxCtgClones * sizeof(struct vertex));
   if(local_a==NULL)
      messcrash("ERROR -- Vertex-number range too large - failed memory alloc\n");

   start=-1;
   sub_info->n=0;
   last=first;

   /* Initialize vertices and set up initial vertex stack, leaving out first*/
   for(i=0; i<MaxCtgClones; i++){
      if((i!=first) && (SMctg[ctgID][i].count==0))
         {local_a[i].count=start; start=i;}
      else
         local_a[i].count=SMctg[ctgID][i].count;
      local_a[i].l=SMctg[ctgID][i].l;
      local_a[i].r=SMctg[ctgID][i].r;
      local_a[i].link=SMctg[ctgID][i].link;
   }

   /* Remove all vertices except those reachable from first*/
   while(start != -1){
      j = start; start = local_a[j].count; p = local_a[j].link;
      while(p != NULL){
         k=p->num;
         local_a[k].count--;
         if(local_a[k].count==0){
            local_a[k].count = start; start = k;
         }
         p=p->next;
      }
   }

   /*Visit remaining vertices and retrieve the one with greatest right end value*/
   local_a[first].count = -1;
   start = first;
   while(start != -1){
      j=start; start = local_a[j].count; p = local_a[j].link;

      /* Visit */
      if(local_a[j].r > max_right){
         sub_info->r_clone= j;
         max_right = local_a[j].r;
      }
      sub_info->clones[sub_info->n]=j;
      sub_info->n++;

      while(p != NULL){
         k=p->num;
         local_a[k].count--;
         if(local_a[k].count==0){
            local_a[k].count = start; start = k;
         }
         p=p->next;
      }
   }
   free(local_a);
   local_a = 0;
}

/********************************
 Allocates memory for the subgraph info struct as needed.
*********************************/
void get_subgraph_mem(struct subgraph **s, int n, int *max){

   if(n >= (*max)-1){
      if(*max==0){
         *max=50;
         *s = (struct subgraph*) calloc(sizeof(struct subgraph),*max);
         if(*s==NULL)  printf("ERROR -- requested too much memory\n");
      }
      else{
         *max+=25;
         *s = (struct subgraph*) realloc(*s, sizeof(struct subgraph)*(*max));
         if(*s==NULL)  printf("ERROR -- requested too much memory\n");
      }
   }
}

/********************************
For a given ctg, find the beginnings and ends of all expressways.
Return number of expressways via num_subgraphs,
*********************************/
void
get_subgraphs(int ctgID, struct subgraph **s, int *num_subgraphs)
{
   int i;
   int max_subgraphs=0;

   for (i = 0; i < num_clones[ctgID]/*MaxCtgClones*/; i++) {
       CLONE *cln = arrp(acedata, clonekey[ctgID][i], CLONE);
       // SD clones don't overlap with anything, but we also shouldn't
       // waste resources by making them into subgraphs
       if (cln->seqstat != SD) {
           if (SMctg[ctgID][i].count == 0) {
               if (*num_subgraphs >= max_subgraphs - 1)/*If no memory, alloc.*/
                   get_subgraph_mem(s, *num_subgraphs, &max_subgraphs);
               (*s)[*num_subgraphs].l_clone = i;
               scan_subgraph(ctgID, i, &((*s)[*num_subgraphs]));
               (*num_subgraphs)++;
           }
       }
   }
}

void
get_mandatory(int ctgID, GPtrArray *mandatory)
{
    struct subgraph *sbgrph = NULL;
    GPtrArray *clones; // holds pointers to FPC CLONEs
    guint i;
    CLONE *prev_cln = NULL;


    clones = g_ptr_array_new();
    for (i = 0; i < num_clones[ctgID]; ++i) {
        CLONE *cln = arrp(acedata, clonekey[ctgID][i], CLONE);
        if (use_mandatory(cln)) {
            g_ptr_array_add(clones, cln);
        }
    }

    // Sort clones, and assemble into subgraphs.
    g_ptr_array_sort(clones, (GCompareFunc)compare_left);
    for (i = 0; i < clones->len; ++i) {
        CLONE *cln = g_ptr_array_index(clones, i);
        int id = get_clone_id(ctgID, cln->clone);
        struct pair *pair = NULL;
        if (prev_cln != NULL)
            pair = find_pair(ctgkey[ctgID][0],
                                       prev_cln->clone, cln->clone);
        if (pair == NULL) {
            sbgrph = g_new0(struct subgraph, 1);
            sbgrph->n = 0;
            sbgrph->l_clone = id;
            sbgrph->r_clone = id;
            g_ptr_array_add(mandatory, sbgrph);
        }
        g_assert(sbgrph != NULL);
        sbgrph->clones[sbgrph->n++] = id;
        g_assert(sbgrph->n <= MAXCTGCLONES);
        sbgrph->r_clone = id;
        prev_cln = cln;
    }

    g_ptr_array_free(clones, TRUE);
    clones = 0;
}

/********************************
 Shortest path function.  Selects vertex u that minimizes dist[u] over
 all u such that inTheTree[u]=false.
 Note: Make sparse matrix a priority queue to correspond to complexity
       given in GR paper.
*********************************/
int
get_next_vertex(EdgeWeight dist[], int inTheTree[], int n)
{
    int u = -1, v;
    EdgeWeight max;
    EdgeWeight* min = &max;

    weight_make_max(&max);
    for (v = 0; v < n; v++) {
        if (!inTheTree[v] && weight_compare(&dist[v], min) < 0) {
            u = v;
            min = &dist[v];
        }
    }
    return u;
}

int get_size_penalty(int sizeL, int sizeR){
  int min, max;

  min=Proj.avginsertsize * 1.3;
  max=Proj.avginsertsize * 2;

  if(min == 0 && max == 0){
    min=MINPAIRSIZE;
    max=MAXPAIRSIZE;
  }

  return CLAMP(-((float)max_weight/(float)(max-min)) *
               (float)(sizeR+sizeL - min) + max_weight,
               0, max_weight);
}



/********************************
 Shortest path function. Computes weight based on fpcOlap.
 Used in weighting edges between expressways.
*********************************/
int get_mtp_weight(int fpcOlap, int fpcSpan){
   float k;   /*konstant*/
   int weight_olap, weight_span, weight_total;

   /*Compute overlap weight*/
   if(fpcOlap<=BESTOLAP){
      k=(float)MAXWEIGHT_OLAP/(float)BESTOLAP;
      weight_olap= -k*fpcOlap + BESTOLAP*k;
   }
    else{
      k=(float)MAXWEIGHT_OLAP/((float)MAXOLAP-(float)BESTOLAP);
      weight_olap=k*fpcOlap - BESTOLAP*k;
      if(weight_olap>MAXWEIGHT_OLAP)
         weight_olap=MAXWEIGHT_OLAP;
   }

   /*Compute span weight*/
   if(fpcSpan<BESTSPAN){
      k=(float)MAXWEIGHT_SPAN/(float)BESTSPAN;
      weight_span= -k*fpcSpan + BESTSPAN*k;
   }
   else{
      weight_span=0;
   }

   weight_total = P1 * weight_olap + P2 * weight_span;
   return weight_total;
}

/********************************
 Shortest path function.  Finds clone c in clones
*********************************/
int get_sp_clone_id(int c, int n, int clones[]){
   int i;

   for(i=0;i<n;i++){
      if(c==clones[i])
         return i;
   }
   return -1;
}

/********************************
 Shortest path function. Runs Dijkstra's shortest path algorithm.
*********************************/


void
dijkstra(int ctgID, int l, int n, int clones[], int parent[],
         EdgeWeight dist[])
{
    int u, v, stage;
    EdgeWeight w;
    int *inTheTree;
    struct edge *edges;
    gboolean prefer_large_p = (prefer_large ? TRUE : FALSE);
    int nedge;

    inTheTree = (int *)malloc(n * sizeof(int));
    if (inTheTree == NULL)
        messcrash("ERROR -- Vertex-number range too large\n");

    for (v = 0; v < n; v++) {
        weight_make_max(&dist[v]);
        inTheTree[v] = 0;
        parent[v] = -1;
    }

    weight_make_min(&dist[l]);
    for (stage = 0; stage < n; stage++) {
        u = get_next_vertex(dist, inTheTree, n);
        if (u != -1) {
            guint u_sm;
            inTheTree[u] = 1;
            u_sm = clones[u];
            edges = SMctg[ctgID][u_sm].link;

            nedge = 0;
            while (edges != NULL) {
                v = get_sp_clone_id(edges->num, n, clones);

                if (v != -1 && !inTheTree[v]) {
#ifdef DIJK_DEBUG
                printf("DIJK:check edge %s %s\n",
                  arrp(acedata, clonekey[ctgID][edges->num], CLONE)->clone,  
                  arrp(acedata, clonekey[ctgID][u_sm], CLONE)->clone); fflush(0);
#endif
                    weight_init(&w, edges->overlap_score, edges->pair_type,
                                prefer_large_p, edges->sizeL, edges->sizeR);
                    nedge++;
                    weight_add(&w, &dist[u]);

                    if (weight_compare(&w, &dist[v]) < 0) { 
                        weight_copy(&dist[v], &w);
                        parent[v] = u;
#ifdef DIJK_DEBUG
                        printf("DIJK:parent of %s is %s with weights bss:%d fp:%d\n",
                          arrp(acedata, clonekey[ctgID][edges->num], CLONE)->clone,  
                          arrp(acedata, clonekey[ctgID][u_sm], CLONE)->clone,  
                            (int)w.bss,(int)w.fp); fflush(0);
#endif
                    }
                    else if (parent[v] == -1) {
#ifdef DIJK_DEBUG
                        printf("WARNING: %s no valid links\n",arrp(acedata, clonekey[ctgID][edges->num], CLONE)->clone);
                        fflush(0);
#endif 
                    }
                }
                edges = edges->next;
            }
#ifdef DIJK_DEBUG
            if (nedge == 0) {
                printf("DIJK:%s has no edges\n",
                  arrp(acedata, clonekey[ctgID][u_sm], CLONE)->clone); fflush(0);
            }
#endif
        }
    }

    free(inTheTree);
    inTheTree = 0;
}

/********************************
 Results function.  Copies next available expressway name into name.
*********************************/
void get_exp_name(char *name){
   sprintf(name,"exp%d",exp_num);
   exp_num++;
}

/********************************
 Results function.  Adds results for current expressway/mtp to results struct
*********************************/
static void
add_to_results(int ctgID, struct subgraph *sgf,
               int *parent, guint left, guint right)
{
    struct expdata expd;
    struct cloneinfo *clonePtr;
    struct cloneinfo tempPtr;
    char c[CLONE_SZ];
    int i, j, max_len, len;
    int num_sp_clones;
    struct ctgdata **results;
    int p_index, n_index, c_index, index;
    char tmp[CLONE_SZ];
    CLONE *clp, *max_clp;
    struct contig *cptr;


    results = exp_results;

    /*Only happens once for each contig*/
    if (results[ctgID]==NULL) {
        /* Contig info */
        results[ctgID] = (struct ctgdata *)malloc(sizeof(struct ctgdata));
        if (results[ctgID] == NULL)
            messcrash("ERROR -- too many ctgs\n");
        results[ctgID]->ctg = ctgkey[ctgID][0];
        results[ctgID]->start = ctgkey[ctgID][1];
        results[ctgID]->end = ctgkey[ctgID][2];
        results[ctgID]->expressways = g_array_new(FALSE, FALSE,
                                                  sizeof(struct expdata));
    }

    get_exp_name(expd.name);
    if(sgf != NULL) {
    expd.clones =
        (struct cloneinfo *)malloc(sgf->n * sizeof(struct cloneinfo));
    if (expd.clones == NULL)
        messcrash("ERROR -- too many clones MTP\n");
    clonePtr = expd.clones;
    expd.l = SMctg[ctgID][left].l;
    expd.r = SMctg[ctgID][right].r;

    /* Write clones in shortest path, from last to first */
    index = clonekey[ctgID][right];
    strcpy(c, array(acedata, index, CLONE).clone);
    g_assert(arrp(acedata, index, CLONE)->seqstat != SD);

    index = clonekey[ctgID][left];
    strcpy(tmp, array(acedata, index, CLONE).clone);
    g_assert(arrp(acedata, index, CLONE)->seqstat != SD);

    p_index = get_sp_clone_id(right, sgf->n, sgf->clones);
    for (i = 0; strcmp(c, tmp); i++) {
        int cidx;
        g_assert(i < sgf->n);
        strcpy(clonePtr[i].name, c);
        clonePtr[i].l = SMctg[ctgID][get_clone_id(ctgID, c)].l;
        clonePtr[i].r = SMctg[ctgID][get_clone_id(ctgID, c)].r;
        cidx = get_sp_clone_id(get_clone_id(ctgID, c),
                               sgf->n, sgf->clones);
        g_assert(cidx == p_index);
        p_index = parent[cidx];
        g_assert(p_index != -1);
        n_index = sgf->clones[p_index];
        c_index = clonekey[ctgID][n_index];
        g_assert(arrp(acedata, c_index, CLONE)->seqstat != SD);
        strcpy(c, array(acedata, c_index, CLONE).clone);
    }
    strcpy(clonePtr[i].name, c);   /*Put first clone in array.*/
    clonePtr[i].l = SMctg[ctgID][get_clone_id(ctgID, c)].l;
    clonePtr[i].r = SMctg[ctgID][get_clone_id(ctgID, c)].r;

    num_sp_clones = i + 1;
    expd.n = num_sp_clones;

    /* Reverse the list */
    j = num_sp_clones-1;
    for (i = 0; i < num_sp_clones / 2; i++) {
        strcpy(tempPtr.name, clonePtr[i].name);
        tempPtr.l = clonePtr[i].l;
        tempPtr.r = clonePtr[i].r;
        strcpy(clonePtr[i].name, clonePtr[j].name);
        clonePtr[i].l = clonePtr[j].l;
        clonePtr[i].r = clonePtr[j].r;
        strcpy(clonePtr[j].name, tempPtr.name);
        clonePtr[j].l = tempPtr.l;
        clonePtr[j].r = tempPtr.r;
        j--;
    }
    } else {
      expd.clones =
        (struct cloneinfo *)malloc(sizeof(struct cloneinfo));
      if (expd.clones == NULL)
        messcrash("ERROR -- too many clones MTP\n");
      clonePtr = expd.clones;
      find_Clone(ctgkey[ctgID][0]);
      orderclones(ctgkey[ctgID][0]);
      max_len = 0;
      for(cptr=root; cptr->new!=NULL; cptr=cptr->new) {
        clp=arrp(acedata, cptr->next, CLONE);
        len = clp->y - clp->x + 1;
        if(len > max_len)  {
          max_len = len;
          max_clp = clp;
        }
      }
      expd.l = max_clp->x;
      expd.r = max_clp->y; 
      clonePtr[0].l = max_clp->x;
      clonePtr[0].r = max_clp->y;
      strcpy(clonePtr[0].name, max_clp->clone);
      expd.n = 1;
    }
    g_array_append_val(results[ctgID]->expressways, expd);
}

/********************************
 Sorts the expressways by span.
*********************************/
int sortExp(const struct expdata *x, const struct expdata *y){
   int span_x, span_y;

   span_x=x->r - x->l;
   span_y=y->r - y->l;
   if(span_x < span_y) return 1;
   if(span_x > span_y) return -1;
   return 0;
}

/********************************
 Writes the contents of results to file for reading by pairs2.pl
cari - this appears not to be used, remove?
*********************************/
void
write_results(int ctgs)
{
    int i,j,k;
    char filName[256];
    FILE *fpout;
    struct ctgdata **results;

    results = exp_results;

    for (i = 0; i < ctgs; i++) {
        if (results[i] == NULL)
            continue;   /*Not all ctgs have expressways*/
        g_array_sort(results[i]->expressways, (GCompareFunc)sortExp);
    }
    sprintf(filName, "MTPdata/pick.exp");
    fflush(stdout);
    printf("Writing %s...", filName);
    if ((fpout = fopen(filName, "w")) == NULL)
        printf("Error opening %s\n", filName);
    for (i = 0; i < ctgs; i++) {
        if (results[i] == NULL)
            continue;   /*Not all ctgs have expressways*/
        fprintf(fpout,"ctg%d span: %d to %d\n", results[i]->ctg,
                results[i]->start, results[i]->end);
        for (j = 0; j < results[i]->expressways->len; j++) {
            struct expdata *expd;
            expd = &g_array_index(results[i]->expressways, struct expdata, j);
            fprintf(fpout,
                    "%s Span: %d to %d Positive olap: 0/0/0 "
                    "Negative olap: 0/0/0\n",
                    expd->name, expd->l, expd->r);
            for (k = 0; k < expd->n; k++) {
                fprintf(fpout, "%s %d %d\n", expd->clones[k].name,
                        expd->clones[k].l, expd->clones[k].r);
            }
            fprintf(fpout, "\n");
        }
    }
    fclose(fpout);
    printf("done.\n");
}

/********************************
 Write an fpc file with the expressways as clones with remarks saying which
 clones the expressways contain.
*********************************/
void
write_fpc(int ctgs)
{
    int i, j, k;
    char filName[256];
    FILE *fpout;
    int clones = 0;
    struct ctgdata **results;

    results = exp_results;
    sprintf(filName, "MTPdata/expr.fpc");
    printf("Writing %s...", filName);
    fflush(stdout);
    if ((fpout = fopen(filName,"w")) == NULL)
        printf("Error opening %s\n",filName);
    /* Write project header */
    for (i = 0; i < ctgs; i++) {
        /*Get total number of clones (that is, expressways)*/
        if (results[i] == NULL) {
            /*Next ctg if no expressways*/
            clones++;
            continue;
        }
        clones += results[i]->expressways->len + 1;  /* Add 1 clone for span */
    }
    fprintf(fpout, "// fpc project expr\n");
    fprintf(fpout, "// Contigs %d  Clones %d\n\n", ctgs, clones);

    for (i = 0; i < ctgs; i++){
        /* Make a clone for the ctg span */
        fprintf(fpout, "Clone : \"CTG%d_SPAN\"\n", ctgkey[i][0]);
        fprintf(fpout, "Map \"ctg%d\" Ends Left %d.000\n",
                ctgkey[i][0], ctgkey[i][1]);
        fprintf(fpout, "Map \"ctg%d\" Ends Right %d.000\n\n",
                ctgkey[i][0], ctgkey[i][2]);
        if (results[i] == NULL)
            continue;   /*Next ctg if no expressways*/
        for (j = 0; j < results[i]->expressways->len; j++) {
            struct expdata *expd;
            expd = &g_array_index(results[i]->expressways, struct expdata, j);
            /* Make a clone for each expressway */
            fprintf(fpout, "Clone : \"%s\"\n", expd->name);
            fprintf(fpout, "Map \"ctg%d\" Ends Left %d.000\n",
                    results[i]->ctg, expd->l);
            fprintf(fpout,"Map \"ctg%d\" Ends Right %d.000\n",
                    results[i]->ctg, expd->r);
            for (k = 0; k < expd->n; k++) {
                fprintf(fpout, "Fpc_remark \"%s\"\n",
                        expd->clones[k].name);
            }
            fprintf(fpout, "\n");
        }
    }
    fclose(fpout);
    printf("done.\n");
}

/********************************
 Frees memory that had been alloc'd for a.
*********************************/
void free_a(int ctgs){
   int i,j;
   struct edge *p1, *p2;

   for(i=0;i<ctgs;i++){
    if (SMctg[i] != NULL) {
      for(j=0;j<num_clones[i];j++){
         p1=SMctg[i][j].link;
         while(p1!=NULL){
            p2=p1->next;
            free(p1);
            p1=p2;
         }
      }
      free(SMctg[i]);
      SMctg[i] = 0;
     }
   }
}

/********************************
 Frees memory that had been alloc'd for mtp.
*********************************/
void free_mtp(int ctgs){
  int i;

  for(i=0;i<ctgs;i++){
    if(mtp[i]==NULL) continue;
    if (mtp[i]->covered != NULL) {
        free(mtp[i]->covered);
        mtp[i]->covered = 0;
    }
    free(mtp[i]);
    mtp[i] = 0;
  }
}

/********************************
 Frees memory that had been alloc'd for results.
*********************************/
void
free_results(int ctgs)
{
   int i, j;

   for (i = 0; i < ctgs; i++) {
      if (exp_results[i] == NULL)
          continue;
      for (j = 0; j < exp_results[i]->expressways->len; j++) {
          struct expdata *expd;
          expd = &g_array_index(exp_results[i]->expressways, struct expdata,
                                j);
          if (expd->clones != NULL) {
              free(expd->clones);
              expd->clones = 0;
          }
      }
      g_array_free(exp_results[i]->expressways, TRUE);
      exp_results[i]->expressways = 0;
      free(exp_results[i]);
      exp_results[i] = 0;
   }
}

/********************************
 Sorts clones in MTP by left values of expressways, followed by
 left values of clones, in ascending order.
*********************************/
int sortMtp(struct pickedcloneinfo *x, struct pickedcloneinfo *y){
   if(x->exp->l > y->exp->l) return 1;
   if(x->exp->l < y->exp->l) return -1;
   if(x->l > y->l) return 1;
   if(x->l < y->l) return -1;
   return 0;
}

/********************************
 Returns the percent of the current expressway already covered by existing expressways.
*********************************/
int get_percent_olap(int *covered, int n, int l, int r){
   int i;
   int count=0;
   int pct=0;

   if(l<0)
       printf("WARNING: left end < 0\n");
   if(r>n)
       printf("WARNING: right end > n (n=%d, r=%d)\n", n, r);
   for (i = MAX(0, l); i <= MIN(n, r); i++){
      if(covered[i]) count++;
   }
   pct = (float)((float)count/(float)(r-l+1)) * 100;
   return pct;
}

/********************************
 Fill mtp structure with clones greedily selected to be in the MTP.
*********************************/
void
select_mtp_clones(int ctgs, GString *warn_string)
{
    int i, j, k, p;
    struct expdata *expPtr;
    int ctg_span;
    int ctgnum;

    for (i = 0; i < ctgs; i++) {
        struct mtpdata *mtpd;

        ctgnum = ctgkey[i][0];

        if (exp_results[i] == NULL) {
            mtp[i] = NULL;
            continue;
        }
        mtp[i] = (struct mtpdata *)malloc(sizeof(struct mtpdata));       
        if (mtp[i] == NULL)
            messcrash("ERROR -- too many ctgs\n");

        /* WMN 10/10/06 dynamically allocate this to reduce memory */
        mtp[i]->picked_clones = (struct pickedcloneinfo*)malloc((contigs[ctgnum].count + 1)*sizeof(struct pickedcloneinfo));

        mtpd = mtp[i];
        ctg_span = exp_results[i]->end - exp_results[i]->start + 1;

        mtpd->covered = (int *)malloc(sizeof(int) * (ctg_span + 1));
        if (mtpd->covered == NULL)
            messcrash("ERROR -- contig too large\n");
        for (j = 0; j < ctg_span; j++)
            mtpd->covered[j] = 0;
        mtpd->span = ctg_span;
        mtpd->n = 0;

        /* Greedily cover contigs with expressways.*/
        for (k = 0; k < exp_results[i]->expressways->len; k++) {
            expPtr = &g_array_index(exp_results[i]->expressways,
                                    struct expdata, k);
            /* Get clones for this expressway*/
            for (p = 0; p < expPtr->n; p++) {
                int cln_index;
                struct pickedcloneinfo *picked;
                CLONE *cln;
                g_assert(mtpd->n < MAXCTGCLONES);
                picked = &(mtpd->picked_clones[mtpd->n]);
                if (!fppFind(acedata, expPtr->clones[p].name, &cln_index,
                             cloneOrder)) {
                    g_error("FPC ERROR could not find clone %s \n",
                            expPtr->clones[p].name);
                }
                picked->l = expPtr->clones[p].l;
                picked->r = expPtr->clones[p].r;
                picked->exp = expPtr;
                picked->pos_in_exp = p + 1;
                picked->total_pos = expPtr->n;
                picked->excess = 0;

                picked->index = cln_index;
                cln = arrp(acedata, cln_index, CLONE);
                picked->mandatory = ( use_mandatory(cln) ? TRUE : FALSE);
                mtpd->n++;
                g_assert(mtpd->n <= contigs[ctgnum].count); /* WMN 10/10/06 */
            }
            for (p = expPtr->l - ctgkey[i][1];
                 p <= expPtr->r - ctgkey[i][1]; p++) {
                g_assert(p <= ctg_span);
                mtpd->covered[p] = 1;
            }
        }
        qsort(mtpd->picked_clones, mtpd->n,
              sizeof(struct pickedcloneinfo),(void *)sortMtp);

        for (p = 1; p < mtpd->n - 1; ++p) {
            if (mtpd->picked_clones[p - 1].mandatory &&
                !mtpd->picked_clones[p].mandatory &&
                mtpd->picked_clones[p + 1].mandatory) {
                g_string_sprintfa(warn_string,
                                  "Single clone (%s) in contig %d "
                                  "bridges two sequenced or tile clones "
                                  "(%s and %s).\n",
                                  clone_idx2name(mtpd->picked_clones[p].index), ctgkey[i][0],
                                  clone_idx2name(mtpd->picked_clones[p - 1].index),
                                  clone_idx2name(mtpd->picked_clones[p + 1].index));
            }
        }
    }
}

/********************************
 Marks clones that can be removed without altering the contiguity of a mtp.
*********************************/
void
mark_excess(int ctgs)
{
    int i, j, k;
    struct pickedcloneinfo *x, *y, *z;
    int rem_x, rem_y;  /*Number of clones removable from x, y respectively*/
    int trimmed_x_olap, trimmed_y_olap;  /*Overlap of x, y after
                                          * excess clones removed*/
    int weight_x, weight_y;  /*Weight of overlap of x, y after excess
                              * clones removed*/

    for (i = 0; i < ctgs; i++) {
        if (mtp[i] == NULL)
            continue;
        for (j = 1; j < mtp[i]->n; j++) {
            x = &mtp[i]->picked_clones[j-1];
            y = &mtp[i]->picked_clones[j];
            if (strcmp(x->exp->name, y->exp->name)){   /*Junction point*/
                if (x->r - MINOLAP > y->l){  /*We have overlap*/
                    /* Try removing from x*/
                    k = 0;
                    z = x;
                    while (!z->mandatory && (z->r - MINOLAP > y->l) &&
                           (j - k - 1 > 0)) {
                        // current z is removable, advance to next
                        z = &mtp[i]->picked_clones[j - (++k) - 1];
                    }
                    trimmed_x_olap = z->r - y->l + 1;
                    rem_x = k;

                    /* Try removing from y*/
                    k = 0;
                    z = y;
                    while (!z->mandatory && (z->l + MINOLAP < x->r) &&
                           (j + k + 1 < mtp[i]->n)) {
                        // current z is removable, advance to next
                        z = &mtp[i]->picked_clones[j + (++k)];
                    }
                    trimmed_y_olap = x->r - z->l + 1;
                    rem_y = k;

                    /* Get weight of overlap just like in computing sp's.*/
                    weight_x = get_mtp_weight(trimmed_x_olap, BESTSPAN);
                    weight_y = get_mtp_weight(trimmed_y_olap, BESTSPAN);

                    if (weight_x <= weight_y){   /*Remove from x*/
                        for (k = 0; k < rem_x; k++) {
                            mtp[i]->picked_clones[j-k-1].excess = 1;
                        }
                    }
                    else {
                        for (k = 0; k < rem_y; k++) {
                            mtp[i]->picked_clones[j+k].excess = 1;
                        }
                    }
                }
            }
        }
    }
}

/********************************
 Writes a remark file for all clones picked to be in the MTP.
*********************************/
void write_fpc_remarks(int ctgs){
   int i,j;
   FILE *fpout;
   char filName[256];

   sprintf(filName,"MTPdata/remarks.ace");
   if((fpout=fopen(filName,"w"))==NULL) printf("Error opening %s\n",filName);

   for(i=0;i<ctgs;i++){
      if(mtp[i]==NULL) continue;
      for(j=0;j<mtp[i]->n;j++){
         fprintf(fpout,"\nClone : \"%s\"\n",clone_idx2name(mtp[i]->picked_clones[j].index));
         if(!mtp[i]->picked_clones[j].excess){
            fprintf(fpout,"Remark \"%s_%d/%d\"\n",
                 mtp[i]->picked_clones[j].exp->name,
                 mtp[i]->picked_clones[j].pos_in_exp,
                 mtp[i]->picked_clones[j].total_pos);
         }
         else{
            fprintf(fpout,"Remark \"%s_%d/%dx\"\n",
                 mtp[i]->picked_clones[j].exp->name,
                 mtp[i]->picked_clones[j].pos_in_exp,
                 mtp[i]->picked_clones[j].total_pos);
         }
      }
   }
   fclose(fpout);
}

/********************************
 Writes a file for all clones picked to be in the MTP.
 Added 6/23/03 (fred)
*********************************/
void write_mtp_results(int ctgs){
   int i,j;
   FILE *fpout;
   char filName[256];
   char exp_name[50]="";

   sprintf(filName,"MTPdata/pick.mtp");
   if((fpout=fopen(filName,"w"))==NULL) printf("Error opening %s\n",filName);

   for(i=0;i<ctgs;i++){
      if(exp_results[i]==NULL) continue;
      fprintf(fpout,"ctg%d span: %d to %d\n", exp_results[i]->ctg,
              exp_results[i]->start, exp_results[i]->end);
      if(mtp[i]==NULL) continue;
      for(j=0;j<mtp[i]->n;j++){
        if(strcmp(exp_name, mtp[i]->picked_clones[j].exp->name)){
          /* New expressway*/
          strcpy(exp_name, mtp[i]->picked_clones[j].exp->name);
          fprintf(fpout,"\n%s Span: %d to %d Clones: %d\n", exp_name,
                  mtp[i]->picked_clones[j].exp->l,
                  mtp[i]->picked_clones[j].exp->r,
                  mtp[i]->picked_clones[j].exp->n);
        }
        if(!mtp[i]->picked_clones[j].excess)
	  fprintf(fpout,"%s %d %d\n", clone_idx2name(mtp[i]->picked_clones[j].index),
		  mtp[i]->picked_clones[j].l, mtp[i]->picked_clones[j].r);
        else
	  fprintf(fpout,"%s %d %d *\n", clone_idx2name(mtp[i]->picked_clones[j].index),
		  mtp[i]->picked_clones[j].l, mtp[i]->picked_clones[j].r);
      }
      fprintf(fpout,"\n");
   }
   fclose(fpout);
}

void
load_picks(int ctgs, GString *warn_string)
{
    int i, j;
    int exp_pos1, exp_pos2, exp_cnt, ljunction, rjunction;
    struct pair *picked_pair;
    CLONE *cln1, *cln2;
    long size_sum = 0;
    int size_cnt = 0;

    stats.num_fp_picked = 0;
    stats.num_bss_picked = 0;
    stats.num_mandatory_picked = 0;
    stats.num_exp_jct = 0;

    for (i = 0; i < ctgs; i++) {
        if (exp_results[i] == NULL)
            continue;
        if (mtp[i] == NULL)
            continue;
        ljunction = rjunction = FALSE;
        for (j = 0; j < (mtp[i]->n - 1); j++) {
            if (rjunction){
                rjunction = FALSE;
            }
            exp_pos1 = mtp[i]->picked_clones[j].pos_in_exp;
            exp_pos2 = mtp[i]->picked_clones[j+1].pos_in_exp;
            exp_cnt = mtp[i]->picked_clones[j+1].total_pos;

            if (mtp[i]->picked_clones[j+1].pos_in_exp ==
                mtp[i]->picked_clones[j+1].total_pos)
                rjunction = TRUE;
            else
                rjunction = FALSE;

            if (mtp[i]->picked_clones[j].pos_in_exp == 1)
                ljunction = TRUE;
            else
                ljunction = FALSE;

#ifdef MTP_DEBUG
            printf("PAIR %d %d %d\n", j, mtp[i]->picked_clones[j].l,
                                    mtp[i]->picked_clones[j+1].r);
            fflush(0);
#endif

            picked_pair = find_pair(ctgkey[i][0],
                                    clone_idx2name(mtp[i]->picked_clones[j].index),
                                    clone_idx2name(mtp[i]->picked_clones[j+1].index));
            if (picked_pair == NULL) {

                /* WMN we add pairs for the expressway junctions, to make display consistent */

                picked_pair = (struct pair*)malloc(sizeof(struct pair));
                picked_pair->pair_type = EXP;
                strcpy(picked_pair->c1,clone_idx2name(mtp[i]->picked_clones[j].index));
                strcpy(picked_pair->c2,clone_idx2name(mtp[i]->picked_clones[j+1].index));
                *picked_pair->cspan = 0;
                *picked_pair->clflank = 0;
                *picked_pair->crflank = 0;  
                picked_pair->x1 = mtp[i]->picked_clones[j].l;
                picked_pair->y1 = mtp[i]->picked_clones[j].r,
                picked_pair->x2 = mtp[i]->picked_clones[j+1].l,
                picked_pair->y2 =  mtp[i]->picked_clones[j+1].r,
                picked_pair->index1 = mtp[i]->picked_clones[j].index;   
                picked_pair->index2 = mtp[i]->picked_clones[j+1].index;   
                picked_pair->indexspan = 0;   
                picked_pair->indexlflank = 0;   
                picked_pair->indexrflank = 0; 

                cln1 = arrp(acedata, mtp[i]->picked_clones[j].index, CLONE);
                cln2 = arrp(acedata, mtp[i]->picked_clones[j+1].index, CLONE);
                picked_pair->olap = Proj.avgbandsize*(MIN(cln1->y, cln2->y) - MAX(cln1->x, cln2->x) + 1);
                picked_pair->c1size = (mtp[i]->picked_clones[j].r - mtp[i]->picked_clones[j].l)*Proj.avgbandsize;   
                picked_pair->c1size_from_file = 0;   
                picked_pair->c2size = (mtp[i]->picked_clones[j+1].r - mtp[i]->picked_clones[j+1].l)*Proj.avgbandsize;   
                picked_pair->c2size_from_file = 0; 
                
                printf("Ctg %d: expressway junction between %s and %s\n",ctgkey[i][0],cln1->clone,cln2->clone);fflush(0); 
            }

            add_to_picks(ctgkey[i][0], picked_pair, exp_pos1,
                         exp_pos2, exp_cnt, ljunction, rjunction);


            size_sum += picked_pair->c1size + picked_pair->c2size;
            size_cnt += 2;
            if (picked_pair->pair_type == FINGERPRINT)
                stats.num_fp_picked++;
            else if (picked_pair->pair_type == BSS_DRAFT)
                stats.num_bss_picked++;
            else if (picked_pair->pair_type == MANDATORY)
                stats.num_mandatory_picked++;
            else if (picked_pair->pair_type == EXP)
                stats.num_exp_jct++;
            else
                g_assert_not_reached();
        }
        if(mtp[i]->n == 1) {
          picked_pair = (struct pair*)malloc(sizeof(struct pair));
          picked_pair->pair_type = SNG;
          strcpy(picked_pair->c1,clone_idx2name(mtp[i]->picked_clones[0].index));
          *picked_pair->cspan = 0;
          *picked_pair->clflank = 0;
          *picked_pair->crflank = 0;  
          picked_pair->x1 = mtp[i]->picked_clones[0].l;
          picked_pair->y1 = mtp[i]->picked_clones[0].r,
          picked_pair->x2 = 0;
          picked_pair->y2 = 0; 
          picked_pair->index1 = mtp[i]->picked_clones[0].index;   
          picked_pair->index2 = 0;
          picked_pair->indexspan = 0;   
          picked_pair->indexlflank = 0;   
          picked_pair->indexrflank = 0; 

          picked_pair->olap = 0;
          picked_pair->c1size = (mtp[i]->picked_clones[0].r - mtp[i]->picked_clones[0].l)*Proj.avgbandsize;   
          picked_pair->c1size_from_file = 0;   
          picked_pair->c2size = 0; 
          picked_pair->c2size_from_file = 0; 
          add_to_picks(ctgkey[i][0], picked_pair, 0,
                         0, 1, 0, 0);
          size_sum += mtp[i]->picked_clones[0].r - mtp[i]->picked_clones[0].l + 1;
          size_cnt++;
        }
    }

    if (size_cnt != 0)
        printf("Average MTP clone size: %ld\n", size_sum/size_cnt);
}

/********************************
 Given a ctg and two clone id's, returns their bssOlap
*********************************/
long
get_bssolap(int ctg, int c1, int c2)
{
   struct edge *ptr;
   long result = 0;

   ptr = SMctg[ctg][c1].link;
   while (ptr != NULL && ptr->num != c2)
       ptr = ptr->next;
   if (ptr != NULL && ptr->num == c2)
       result = ptr->bssOlap;
   else {
       result = (MIN(SMctg[ctg][c1].r, SMctg[ctg][c2].r) -
                 MAX(SMctg[ctg][c1].l, SMctg[ctg][c2].l) + 1) *
           Proj.avgbandsize;
   }
   if (result > MAXBASEOLAP) {
      result = MAXBASEOLAP;
   }
   else if (result < -MAXNEGBASEOLAP) {
      result = -MAXNEGBASEOLAP;
   } 

   return result;
}

/********************************
 Fills distribution arrays, averages, etc.
*********************************/
void compute_stats(int ctgs){
   int i,j;
   int p,q;
   struct pickedcloneinfo *x, *y;
   int count=0;
   int left_exp_end, right_exp_end;
   int total_covered=0;
   long bssolap;
   int ctg_gap_total;
   int ctg_gap_left, ctg_gap_right;
   gboolean has_coverage;
   struct mtpdata *mtpd;

   for(i=0;i<MAXBASEOLAP/REPORTBIN;i++){
      stats.exp_pos_cln_olap[i]=0;
      stats.exp_neg_cln_olap[i]=0;
   }
   stats.total_gap=0;
   stats.total_exp_kb=0;
   stats.total_exp_cln=0;
   stats.total_pos_cln_olap=0;
   stats.total_neg_cln_olap=0;
   for(i=0;i<ctgs;i++){
      if(mtp[i]==NULL) {
          stats.ctgstats[i] = NULL;
          continue;
      }
      mtpd = mtp[i];
      ctg_gap_total=0;
      ctg_gap_left = -1;
      ctg_gap_right = 0;
      stats.ctgstats[i] =
         (struct ctgstatdata *)malloc(sizeof(struct ctgstatdata));
      if(stats.ctgstats[i]==NULL)
          messcrash("ERROR -- couldn't alloc memory for stats.\n");
      stats.ctgstats[i]->num_gaps = 0;
      stats.ctgstats[i]->ctg = ctgkey[i][0];
      stats.ctgstats[i]->ctg_len = ctgkey[i][2] - ctgkey[i][1] + 1;
      stats.ctgstats[i]->overlap = 0;

      /* Distribution of gaps */
      count = 0;
      if(mtp[i]->covered[0]) {
          total_covered++;   /*Get first cbunit if covered*/
          ctg_gap_left = 0;
      }
      for(j=1;j<mtp[i]->span;j++){
         p=mtp[i]->covered[j-1];
         q=mtp[i]->covered[j];
         count++;
         if((p!=q) && (!p)){
            if (ctg_gap_left < 0)
                ctg_gap_left = count;
            else {
                stats.ctgstats[i]->num_gaps++;
                stats.total_gap+=(count*Proj.avgbandsize)/1000;
                ctg_gap_total+=count;
            }
            count=0;
         }
         else if(p!=q){
            count=0;
         }
         if (q) {
             total_covered++;
         }
      }
      if (!q)
          ctg_gap_right = count;

      /* Distribution of junction overlaps and picked expressway lengths*/
      count=0;
      has_coverage = FALSE;
      left_exp_end=mtp[i]->picked_clones[0].l;
      for(j=1;j<mtp[i]->n;j++){
         x=&mtp[i]->picked_clones[j-1];
         y=&mtp[i]->picked_clones[j];
         count++;
            /* Picked expressway lengths */
         if((strcmp(x->exp->name, y->exp->name)) || (j==mtp[i]->n-1)){
            has_coverage = TRUE;
            if(j==mtp[i]->n-1){  /* Take care of last clone in last expressway in contig*/
               right_exp_end=y->r;
               stats.total_exp_cln+=count+1;
            }
            else{
               right_exp_end=x->r;
               stats.total_exp_cln+=count;
            }
            stats.total_exp_kb+=((right_exp_end-left_exp_end+1)*Proj.avgbandsize)/1000;
            count=0;
            left_exp_end=y->l;
         }

         bssolap=get_bssolap(i,get_clone_id(i,clone_idx2name(x->index)),get_clone_id(i,clone_idx2name(y->index)));
         if(bssolap>=0){
             g_assert(bssolap/REPORTBIN < G_N_ELEMENTS(stats.exp_pos_cln_olap));
             stats.exp_pos_cln_olap[bssolap/REPORTBIN]++;
             stats.total_pos_cln_olap+=bssolap;
         }
         else{
             g_assert(abs(bssolap)/REPORTBIN <
                      G_N_ELEMENTS(stats.exp_pos_cln_olap));
             stats.exp_neg_cln_olap[abs(bssolap)/REPORTBIN]++;
             stats.total_neg_cln_olap+=abs(bssolap);
         }
         stats.ctgstats[i]->overlap += MAX(bssolap, 0);

      }

      if (has_coverage)
          stats.ctgstats[i]->percent_covered =
              ((float)(stats.ctgstats[i]->ctg_len - ctg_gap_total -
                       ctg_gap_left - ctg_gap_right) /
               (float)stats.ctgstats[i]->ctg_len) * 100;
      else {
          stats.total_exp_cln++;
          stats.ctgstats[i]->percent_covered = 
            ((float)(mtpd->picked_clones[0].r - mtpd->picked_clones[0].l) /
            (float)mtpd->span) *100; 
      }
      stats.ctgstats[i]->gap_len = ctg_gap_total;
   }
   if (total_map_size != 0)
     stats.percent_covered=(float)total_covered/(float)total_map_size * 100;
   else
     stats.percent_covered=0;
}

/********************************
 Prints the contents of the stats to the screen or to a file if outfile != "stdout"
*********************************/
void print_stats(int ctgs, char *outfile){
   FILE *fpout;
   int i,j,k;
   int num;
   struct ctgstatdata *ptr;
   int total_gap=0, total_ctg=0, cnt_mandatory=0, cnt_mtp;
   CLONE *cln;
   char temp[50];
   int pct;
   int total_pairs_picked = 0;
   int clones_picked = 0, single_clones = 0;

   if(!strcmp(outfile,"stdout")) fpout=stdout;
   else if((fpout=fopen(outfile,"w"))==NULL)
        messcrash("Error opening %s\n",outfile);
   fprintf(fpout,"Contig totals: (in CB units)\n");
   fprintf(fpout,"Contig    Ctg len   # MTP  overlap   # of gaps   gap length   %%covered\n");
   fprintf(fpout,"-------   -------   -----  -------   ---------   ----------   --------\n");

   for(i=0;i<ctgs;i++){
      if (num_clones[i] == 0) continue;
      ptr=stats.ctgstats[i];
      if(ptr==NULL || ptr->percent_covered == 0) continue; 

      total_pairs_picked += pick_cnt[ptr->ctg];

      // find number of MTP and show gaps WMN:expressway jncts are added as pairs
        // now so no gaps found
      k = ptr->ctg;
      if(pick_cnt[k] == 1 && picks[k].p[0].exp_cnt == 1) 
        cnt_mtp = 1;
      else {      
        for (temp[0] = '\0', cnt_mtp=1, j=0; j<pick_cnt[k]; j++) {
           if (temp[0]!='\0' &&   // GAP 
               strcmp(temp, picks[k].p[j].picked->c1)!=0) {
               cnt_mtp++; 
               printf("GAP %s\n",temp);
           }
           cnt_mtp++;
           strcpy(temp,  picks[k].p[j].picked->c2);
        } 
      }

       // count mandatory clones used
      for (j = 0; j < mtp[i]->n; j++) {
         cln = arrp(acedata, mtp[i]->picked_clones[j].index, CLONE); 
         if (use_mandatory(cln)) cnt_mandatory++;       
      }
      clones_picked += mtp[i]->n;
      if(mtp[i]->n == 1)
        single_clones++;

      fprintf(fpout,"ctg%-4d   %7d    %4d  %7d   %9d   %10d   %7d%%\n",
           ptr->ctg, ptr->ctg_len, cnt_mtp, 
           UDIV_RND(ptr->overlap, Proj.avgbandsize),
           ptr->num_gaps, ptr->gap_len, ptr->percent_covered);
      total_gap+=ptr->num_gaps;
      total_ctg++;
   }
   fprintf(fpout,"\n");

   fprintf(fpout,"Clone overlap (base pairs):\n");
   fprintf(fpout,"Positive:\n   ");
   num=k=0;
   do{
      for(j=0, i=k;(j<12) && (i<MAXBASEOLAP/REPORTBIN);i++){
         if(stats.exp_pos_cln_olap[i]!=0){
            fprintf(fpout, "%6d-",i*REPORTBIN);
            j++;
         }
      }
      fprintf(fpout,"\n   ");
      for(j=0, i=k;(j<12) && (i<MAXBASEOLAP/REPORTBIN);i++){
         if(stats.exp_pos_cln_olap[i]!=0){
            fprintf(fpout, "%6d ",i*REPORTBIN+REPORTBIN-1);
            j++;
         }
      }
      fprintf(fpout,"\n   ");
      for(j=0, i=k;(j<12) && (i<MAXBASEOLAP/REPORTBIN);i++){
         if(stats.exp_pos_cln_olap[i]!=0){
            fprintf(fpout,"%6d ",stats.exp_pos_cln_olap[i]);
            num+=stats.exp_pos_cln_olap[i];
            j++;
         }
      }
      if(i<MAXBASEOLAP/REPORTBIN)
         fprintf(fpout,"\n\n   ");
      k=i;
   }while(i<MAXBASEOLAP/REPORTBIN);
   fprintf(fpout,"\n   Total positive overlap: %ld\n",
         stats.total_pos_cln_olap);
   if(num!=0) fprintf(fpout,
        "   Average positive overlap: %ld\n",stats.total_pos_cln_olap/num);
   fprintf(fpout,"   Number of positive clone overlaps: %d\n\n",num);

   if (stats.total_neg_cln_olap!=0) {
      fprintf(fpout,"Negative:\n   ");
      num=k=0;
      do{
        j=0;
        for(i=k;(j<6) && (i<MAXNEGBASEOLAP/REPORTBIN);i++){
          if(stats.exp_neg_cln_olap[i]!=0){
             sprintf(temp,"%3d-%-3d ",i*REPORTBIN,i*REPORTBIN+REPORTBIN-1); j++;
             fprintf(fpout,"%-16s",temp);
          }
        }
        if (i == MAXNEGBASEOLAP/REPORTBIN && j < 6) {
         if(stats.exp_neg_cln_olap[i]!=0){
            sprintf(temp,">%3d",i*REPORTBIN); j++;
            fprintf(fpout,"%-8s",temp);
         }       
        }
        fprintf(fpout,"\n        ");
        j=0;
        for(i=k;(j<6) && (i<=MAXNEGBASEOLAP/REPORTBIN);i++){
           if(stats.exp_neg_cln_olap[i]!=0){ fprintf(fpout,"%-16d",stats.exp_neg_cln_olap[i]); 
                   num+=stats.exp_neg_cln_olap[i];j++;}
        }
        if (i == MAXNEGBASEOLAP/REPORTBIN && j < 6) {
           if(stats.exp_neg_cln_olap[i]!=0){
              fprintf(fpout,"%-8d",stats.exp_neg_cln_olap[i]);
              num+=stats.exp_neg_cln_olap[i];
              j++;
           }       
        }

        if(i<MAXNEGBASEOLAP/REPORTBIN)
             fprintf(fpout,"\n\n   ");
        k=i;
      }while(i <= MAXNEGBASEOLAP/REPORTBIN);

      fprintf(fpout,"\n   Total negative overlap: %ld\n",
           stats.total_neg_cln_olap);
      if(num!=0)
         fprintf(fpout,"   Average negative overlap: %ld\n",stats.total_neg_cln_olap/num);
      fprintf(fpout,"   Number of negative clone overlaps (spanned by draft): %d\n\n",num);
   } // end negative overlap

    
   fprintf(fpout, "Clones picked:%d\n",clones_picked);
   if (total_pairs_picked > 0) {
       pct = 100*stats.num_bss_picked/total_pairs_picked;
       fprintf(fpout,"BSS pairs: %d (%d%%)\n", stats.num_bss_picked, pct);
       pct = 100*stats.num_fp_picked/total_pairs_picked;
       fprintf(fpout,"Fingerprint pairs: %d (%d%%)\n", stats.num_fp_picked, pct);
       fprintf(fpout,"Single MTP clones: %d \n", single_clones);
       pct = 100*stats.num_mandatory_picked/total_pairs_picked;
       fprintf(fpout, "Mandatory clone pairs: %d (%d%%)\n", stats.num_mandatory_picked, pct);
       pct = 100*stats.num_exp_jct/total_pairs_picked;
       fprintf(fpout, "Expressway junctions: %d (%d%%)\n", stats.num_exp_jct, pct);


       fprintf(fpout,"\nNumber of clones in MTP: %d\n", stats.total_exp_cln);
       fprintf(fpout, "Number of mandatory clones: %d\n\n", cnt_mandatory);
       fprintf(fpout,"Total gap span: %d kb\n",stats.total_gap);
       fprintf(fpout,"Total MTP span: %d kb\n", stats.total_exp_kb);
       fprintf(fpout,"Percent of map covered: %d%%\n",stats.percent_covered);
   }
   if(strcmp(outfile,"stdout")) fclose(fpout);
}

/********************************
 Writes a text file in a format for entering into an excel spreadsheet.
*********************************/
void write_excel(int ctgs, char *outfile){
   FILE *fpout;
   char filName[256];
   int i,j,k;

   if(!strcmp(outfile,"stdout")) strcpy(filName,"excel.txt");
   else sprintf(filName,"%s.xl.txt",outfile);

   if((fpout=fopen(filName,"w"))==NULL) fprintf(fpout,"Error opening %s\n",filName);

   fprintf(fpout,"Positive clone overlap distribution\n");
   k=0;
   for(i=0;i<MAXBASEOLAP/REPORTBIN;i++){
      if(stats.exp_pos_cln_olap[i]!=0){
         for(j=k+1;j<i;j++) fprintf(fpout,"%d-%d,",j*10,j*10+9);
         fprintf(fpout,"%d-%d,",i*REPORTBIN,i*REPORTBIN+REPORTBIN-1);
         k=i;
      }
   }
   fprintf(fpout,"\n");
   k=0;
   for(i=0;i<MAXBASEOLAP/REPORTBIN;i++){
      if(stats.exp_pos_cln_olap[i]!=0){
         for(j=k+1;j<i;j++) fprintf(fpout,"%d,",stats.exp_pos_cln_olap[j]);
         fprintf(fpout,"%d,",stats.exp_pos_cln_olap[i]);
         k=i;
      }
   }
   fprintf(fpout,"\n\n");

   fprintf(fpout,"Negative clone overlap distribution\n");
   k=0;
   for(i=0;i<MAXBASEOLAP/REPORTBIN;i++){
      if(stats.exp_neg_cln_olap[i]!=0){
         for(j=k+1;j<i;j++) fprintf(fpout,"%d-%d,",j*REPORTBIN,j*REPORTBIN+REPORTBIN-1);
         fprintf(fpout,"%d-%d,",i*REPORTBIN,i*REPORTBIN+REPORTBIN-1);
         k=i;
      }
   }
   fprintf(fpout,"\n");
   k=0;
   for(i=0;i<MAXBASEOLAP/REPORTBIN;i++){
      if(stats.exp_neg_cln_olap[i]!=0){
         for(j=k+1;j<i;j++) fprintf(fpout,"%d,",stats.exp_neg_cln_olap[j]);
         fprintf(fpout,"%d,",stats.exp_neg_cln_olap[i]);
         k=i;
      }
   }
   fprintf(fpout,"\n\n");

   fprintf(fpout,"Contig coverage distribution\n");
   for(i=0;i<ctgs;i++)
      fprintf(fpout,"%d,",ctgkey[i][0]);
   fprintf(fpout,"\n");
   for(i=0;i<ctgs;i++){
      if(stats.ctgstats[i]==NULL)
         fprintf(fpout,"0,");
      else
         fprintf(fpout,"%d,",stats.ctgstats[i]->percent_covered);
   }
   fprintf(fpout,"\n");

   fclose(fpout);
}

/********************************
 Main loop for the expressway picker.
*********************************/
void pick_exp(int output_fpc)
{
    guint i, j, k;
    struct subgraph *subgraphs; /* Contains info for current ctg's subgraphs */
    guint num_subgraphs;
    GArray *left, *right;
    CLONE *cln1;
    CLONE *cln2;

    for (i = 0; i < MAXCTGS; i++) {
        exp_results[i] = NULL;
        /*num_clones[i] = 0;   WMN */
    }

    /*num_ctgs = 0;     WMN */
    left = g_array_new(FALSE, FALSE, sizeof(int));
    right = g_array_new(FALSE, FALSE, sizeof(int));

    load_sparse_matrix();
    for (i = 0; i < num_ctgs; i++) {
        GPtrArray* mandatory;

        GArray* clone_array = g_array_index(clones_ctg, GArray*, ctgkey[i][0]);

        // FIXME: correct boundary offset values
        gint left_bdy = MIN(ctgkey[i][1] + Pz.fromendInt, ctgkey[i][2]);
        gint right_bdy = MAX(ctgkey[i][2] - Pz.fromendInt, ctgkey[i][1]);

        /* WMN 9/11/06 we could probably handle the ==1 case but at the
            moment it hits the asserts below so I excluded it. You
            can't get 1-clone contigs that easily anyway.  */
        if (num_clones[i] == 0 || num_clones[i] == 1)
            continue;
        printf("\rFinding MTP for ctg%d...", ctgkey[i][0]);
        fflush(stdout);
        num_subgraphs = 0;
        subgraphs = NULL;
        //get_subgraphs(i, &subgraphs, &num_subgraphs);

        mandatory = g_ptr_array_new();
        get_mandatory(i, mandatory);

        g_array_set_size(left, 0);
        g_array_set_size(right, 0);

        if (mandatory->len == 0) {
            for (j = 0; j < num_clones[i]; ++j) 
                if (SMctg[i][j].l <= left_bdy && SMctg[i][j].r != ctgkey[i][2]) 
                   g_array_append_val(left, j);
           
            for (j = 0; j < num_clones[i]; ++j) {
                if (SMctg[i][j].r >= right_bdy) {
                  if(SMctg[i][j].r == ctgkey[i][2]) 
                     g_array_append_val(right, j);
                  else {
                   for(k = 0; k < left->len; k++) 
                     if(j == g_array_index(left, int, k)) break;
                   if(k == left->len)
                     g_array_append_val(right, j);
                  }
                }
            }
            /* WMN 9/11/06 for small contigs all clones can wind up left
                or right. In this case just pick the biggest clone.  */
            if (left->len == 0 || right->len == 0) {

                int maxsize = SMctg[i][0].r - SMctg[i][0].l;
                int best_i = 0;
                for (j = 1; j < num_clones[i]; ++j) {
                    if (SMctg[i][j].r - SMctg[i][j].l > maxsize) {
                        best_i = j;
                        maxsize = SMctg[i][j].r - SMctg[i][j].l;
                    }
                }

                add_to_results(i, 0, 0, 0,0);     

            }
            else {
                g_assert(left->len > 0);
                g_assert(right->len > 0);
                span_gap(i, left, right);
            }
        }
        else {
            gint id_left, id_right;
            guint f_idx = 0;
            struct subgraph *sgf;

            // SP from left end of contig to first mandatory subgraph
            sgf = g_ptr_array_index(mandatory, f_idx);
            id_left = sgf->l_clone;
            id_right = sgf->r_clone;

            cln2 = g_array_index(clone_array, CLONE*, id_right);
            cln1 = g_array_index(clone_array, CLONE*, id_left);
#ifdef MTP_DEBUG
            fprintf(stderr,"Span left end to mandatory:left %s right:%s\n",cln1->clone, cln2->clone);
            fflush(0);
#endif

            if (SMctg[i][id_left].l > left_bdy) {
                for (j = 0; j < num_clones[i]; ++j) {
                    if (SMctg[i][j].l <= left_bdy && j != id_left)
                        g_array_append_val(left, j);
                }
                g_assert(left->len > 0);
                g_array_append_val(right, id_left);
                span_gap(i, left, right);
            }
            add_mandatory_to_results(i, sgf);


            // SP's between mandatory subgraphs
            while (++f_idx < mandatory->len) {
                sgf = g_ptr_array_index(mandatory, f_idx);
                id_left = sgf->l_clone;
                g_array_set_size(left, 0);
                g_array_set_size(right, 0);
                g_array_append_val(left, id_right);
                g_array_append_val(right, id_left);

                cln2 = g_array_index(clone_array, CLONE*, id_right);
                cln1 = g_array_index(clone_array, CLONE*, id_left);
#ifdef MTP_DEBUG
                fprintf(stderr,"Span mandatory left:%s right:%s\n",cln1->clone, cln2->clone);
                fflush(0);
#endif

                if ( (cln1->y >= cln2->x && cln1->y <= cln2->y) ||
                        (cln2->y >= cln1->x && cln2->y <= cln1->y) ) {
                    ; // WMN to prevent finding extra mtp between mandatory
                    // clones which overlap in FPC but which aren't part of a pair
                }
                else {
                    span_gap(i, left, right);
                }
                id_right = sgf->r_clone;
                add_mandatory_to_results(i, sgf);
            }

            // SP from last mandatory subgraph to right end of contig
            if (SMctg[i][id_right].r < right_bdy) {
                g_array_set_size(left, 0);
                g_array_set_size(right, 0);
                for (j = 0; j < num_clones[i]; ++j) {
                    if (SMctg[i][j].r >= right_bdy && j != id_right)
                        g_array_append_val(right, j);
                }
                g_assert(right->len > 0);
                g_array_append_val(left, id_right);

                cln2 = g_array_index(clone_array, CLONE*, id_right);
                cln1 = g_array_index(clone_array, CLONE*, id_left);
#ifdef MTP_DEBUG
                fprintf(stderr,"Span mandatory to right end:left %s right:%s\n",cln1->clone, cln2->clone);
                fflush(0);
#endif
                span_gap(i, left, right);
            }
        }

        if (subgraphs!=NULL) {
            free(subgraphs);
            subgraphs = 0;
        }
        g_ptr_array_foreach(mandatory, (GFunc)free_mandatory_subgraph, NULL);
        g_ptr_array_free(mandatory, TRUE);
        mandatory = 0;
    }
    g_array_free(left, TRUE);
    left = 0;
    g_array_free(right, TRUE);
    right = 0;
    printf("\rFinding MTP completed.         \n");
}

/********************************
 Main loop for the mtp picker.
*********************************/
void pick_mtp(char *outfile, int write_xl_file){
   int i;
   GString *warn_string;

   if (num_ctgs > MAXCTGS)
       g_error("Too many contigs!");

   for(i=0;i<MAXCTGS;i++){
      mtp[i]=NULL;
   }
   warn_string = g_string_new("");
   join_expressways(num_ctgs);
   select_mtp_clones(num_ctgs, warn_string);
   mark_excess(num_ctgs);
/*
   write_fpc_remarks(num_ctgs);
   write_mtp_results(num_ctgs);
*/
   load_picks(num_ctgs, warn_string);
   if (warn_string->len > 0) {
       printf("Locations detected on the MTP that may warrant "
              "human inspection:\n%s", warn_string->str);
   }
   compute_stats(num_ctgs);
   print_stats(num_ctgs, outfile); fflush(0);
   if(write_xl_file) write_excel(num_ctgs,outfile);
   /* Clean up */
   g_string_free(warn_string, TRUE);
   warn_string = 0;
   free_results(num_ctgs);
   free_mtp(num_ctgs);
   free_a(num_ctgs);
}

// free_mandatory_subgraph
//
// Free 'struct subgraph' element of a GPtrArray. (Called from
// g_ptr_array_foreach.)
//
static void
free_mandatory_subgraph(struct subgraph *sbgrph, gpointer reserved)
{
    g_free(sbgrph);
    sbgrph = 0;
}

// compare_left
//
// Compare left end coordinates of two clones (qsort-style function).
//
static gint
compare_left(const CLONE **c1, const CLONE **c2)
{
    gint result;

    if ((*c1)->x < (*c2)->x)
        result = -1;
    else if ((*c1)->x > (*c2)->x)
        result = 1;
    else
        result = 0;
    return result;
}

// add_mandatory_to_results
//
// Create an expressway from a subgraph of 'mandatory' clones and add it to
// results.
//
static void
add_mandatory_to_results(guint ctgID, struct subgraph *sgf)
{
    int *sp_parent;
    gint i;
    gap_path_t gp;

    g_assert(sgf->n > 0);
    
#ifdef MTP_DEBUG
    printf("Add mandatory subgraph of size %d to results\n",sgf->n);fflush(0);
#endif
    
    sp_parent = g_new0(gint, sgf->n);
    for (i = 0; i < sgf->n; ++i) {
        sp_parent[i] = i - 1;
#ifdef MTP_DEBUG
        CLONE* clp = arrp(acedata, clonekey[ctgID][sgf->clones[i]], CLONE);
        printf("mandatory:%s\n",clp->clone);fflush(0);
#endif
    }

    gp.sgf = sgf;
    gp.parent = sp_parent;
    gp.sm_left = sgf->clones[0];
    gp.sm_right = sgf->clones[sgf->n - 1];
    gp.sp_left = 0;
    gp.sp_right = sgf->n - 1;
    gp.left_pref = FALSE;
    gp.right_pref = FALSE;
    if (sgf->n > 1) {
        g_assert(path_ok(ctgID, &gp));
    }

    add_to_results(ctgID, sgf, sp_parent, sgf->clones[0],
                   sgf->clones[sgf->n - 1]);
    g_free(sp_parent);
    sp_parent = 0;

}

// span_gap
//
// Create expressways between a set of left and right endpoint clones.
// The idea is to build subgraphs of only the clones inside the gap,
// build expressways from those subgraphs, and then span the gap with
// those expressways.
//
static void
span_gap(guint ctgID, GArray *left_ends, GArray *right_ends)
{
    GArray *gap_subgraphs;
    GArray *gap_clones;
    GArray *spanning_paths;
    guint i;

    gap_clones = find_gap_clones(ctgID, left_ends, right_ends);
    gap_subgraphs =
        find_gap_subgraphs(ctgID, gap_clones, left_ends);
    g_array_free(gap_clones, TRUE);
    gap_clones = 0;
    spanning_paths =
        find_spanning_paths(ctgID, gap_subgraphs, left_ends, right_ends);
    
    trim_and_add_paths(ctgID, spanning_paths);
    
    // a bit of ugliness in gap_path_t exposed here---watch for memory
    // leaks if usage of spanning_paths changes
    for (i = 0; i < spanning_paths->len; ++i) {
        g_free((&g_array_index(spanning_paths, gap_path_t, i))->parent);
        (&g_array_index(spanning_paths, gap_path_t, i))->parent = 0;       
    }
    g_array_free(spanning_paths, TRUE);
    spanning_paths = 0;
    g_array_free(gap_subgraphs, TRUE);
    gap_subgraphs = 0;
}

static GArray*
find_gap_clones(guint ctgID, GArray *left_ends, GArray *right_ends)
{
    struct vertex *sm;
    gint left_bdy, right_bdy;
    guint i;
    GArray *result;

    sm = SMctg[ctgID];

    // Compute left boundary of gap; that is, leftmost left end of
    // clones in 'left_ends'.
    left_bdy = sm[g_array_index(left_ends, gint, 0)].l;
    for (i = 1; i < left_ends->len; ++i) {
        gint l = sm[g_array_index(left_ends, gint, i)].l;
        if (l < left_bdy) {
            left_bdy = l;
        }
    }

    // Compute right boundary of gap; that is, rightmost right end of
    // clones in 'right_ends'.
    right_bdy = sm[g_array_index(right_ends, gint, 0)].r;
    for (i = 1; i < right_ends->len; ++i) {
        gint r = sm[g_array_index(right_ends, gint, i)].r;
        if (r > right_bdy) {
            right_bdy = r;
        }
    }

    // Get all clones inside the gap.
    result = g_array_new(FALSE, FALSE, sizeof(gap_clone_t));
    for (i = 0; i < num_clones[ctgID]; ++i) {
        CLONE *c = arrp(acedata, clonekey[ctgID][i], CLONE);
        gap_clone_t cln;
        if (sm[i].l >= left_bdy && sm[i].r <= right_bdy &&
            c->seqstat != SD) {
            cln.in_subgraph = FALSE;
            cln.sm_index = i;
            cln.left = sm[i].l;
            g_array_append_val(result, cln);
#ifdef MTP_DEBUG
            fprintf(stderr, "GAP CLONE:%s\n", c->clone);
            fflush(0);
#endif
        }
    }
    g_array_sort(result, (GCompareFunc)compare_gap_clones);
    return result;
}

static gint
compare_gap_clones(const gap_clone_t *gc1, const gap_clone_t *gc2)
{
    gint result;

    if (gc1->left < gc2->left)
        result = -1;
    else if (gc1->left > gc2->left)
        result = 1;
    else
        result = 0;
    return result;
}

static GArray*
find_gap_subgraphs(guint ctgID, GArray *gap_clones, GArray *left_ends)
{
    guint i;
    struct vertex *sm;
    GArray *result;

    sm = SMctg[ctgID];

    result = g_array_new(FALSE, FALSE, sizeof(struct subgraph));

    // Get a subgraph for every clone in 'left_ends'.
    for (i = 0; i < left_ends->len; ++i) {
        guint j;
        guint left_sm = g_array_index(left_ends, int, i);
        for (j = 0; j < gap_clones->len; ++j) {
            gap_clone_t *gc = &g_array_index(gap_clones, gap_clone_t, j);
            if (left_sm == gc->sm_index) {
                scan_gap_subgraph(sm, gap_clones, gc, result);
            }
        }
    }

    // Pick up other subgraphs by searching for clones in the gap that
    // are not yet part of any subgraph. Because the list of clones in
    // 'gap_clones' is sorted by the left end of the clones, we're
    // always starting a new subgraph as far to the left in the gap as
    // possible (among those clones not in a subgraph at the time they
    // are checked).
    for (i = 0; i < gap_clones->len; ++i) {
        gap_clone_t *gc = &g_array_index(gap_clones, gap_clone_t, i);
        if (!gc->in_subgraph) {
            scan_gap_subgraph(sm, gap_clones, gc, result);
        }
    }

    g_array_sort_with_data(result, (GCompareDataFunc)compare_subgraphs_left,
                           sm);
    return result;
}

static void
scan_gap_subgraph(const struct vertex *sm, GArray *gap_clones,
                  gap_clone_t *start, GArray *subgraphs)
{
    guint i;
    GPtrArray *subgraph_clones;
    GSList *branches = NULL;

    for (i = 0; i < gap_clones->len; ++i) {
        (&g_array_index(gap_clones, gap_clone_t, i))->in_this_subgraph = FALSE;
    }

    subgraph_clones = g_ptr_array_new();
    start->in_this_subgraph = TRUE;
    branches = g_slist_prepend(branches, start);

    while (branches != NULL) {
        gap_clone_t *branch_clone = branches->data;
        struct edge *edge;
        // Pop head of 'branches', copy it to 'subgraph_clones', and
        // set the clone's 'in_subgraph' field.
        branches = g_slist_remove(branches, branch_clone);
        branch_clone->in_subgraph = TRUE;
        g_ptr_array_add(subgraph_clones, branch_clone);
        // Add to 'branches' all clones connected to 'branch_clone' by
        // an edge.
        edge = sm[branch_clone->sm_index].link;
        while (edge != NULL) {
            gap_clone_t *neighbor = NULL;
            for (i = 0; neighbor == NULL && i < gap_clones->len; ++i) {
                gap_clone_t *gc;
                gc = &g_array_index(gap_clones, gap_clone_t, i);
                if (edge->num == gc->sm_index) {
                    neighbor = gc;
                }
            }
            if (neighbor != NULL && !neighbor->in_this_subgraph) {
                neighbor->in_this_subgraph = TRUE;
                branches = g_slist_prepend(branches, neighbor);
            }
            edge = edge->next;
        }
    }

    // Convert accumulated list of clones into 'struct subgraph' instance.
    if (subgraph_clones->len > 1) { // TODO: >1 or >0?
        struct subgraph sgf;
        guint i;
        sgf.n = 0;
        g_assert(subgraph_clones->len <= MAXCTGCLONES);
        for (i = 0; i < subgraph_clones->len; ++i) {
            gap_clone_t *gc;
            gc = g_ptr_array_index(subgraph_clones, i);
            if (sgf.n == 0) {
                sgf.l_clone = gc->sm_index;
                sgf.r_clone = gc->sm_index;
            }
            else {
                if (sm[gc->sm_index].l < sm[sgf.l_clone].l) {
                    sgf.l_clone = gc->sm_index;
                }
                if (sm[gc->sm_index].r > sm[sgf.r_clone].r) {
                    sgf.r_clone = gc->sm_index;
                }
            }
            sgf.clones[sgf.n++] = gc->sm_index;
        }
        g_array_append_val(subgraphs, sgf);
    }
    g_ptr_array_free(subgraph_clones, TRUE);
}

static GArray*
find_spanning_paths(guint ctgID, GArray *gap_subgraphs,
                    GArray *left_ends, GArray *right_ends)
{
    guint i;
    EdgeWeight *dist = NULL;
    GArray *gap_paths;
    GSList *path_candidates;
    guint path_idx;
    struct vertex *sm;
    GArray *result;
    int cont;

    sm = SMctg[ctgID];
    result = g_array_new(FALSE, FALSE, sizeof(gap_path_t));

    // Compute shortest path going between extremities of each
    // subgraph. The exception is for the graphs containing clones in
    // left_ends or right_ends, for which we prefer to use one of
    // those clones as an endpoint.
    gap_paths = g_array_new(FALSE, FALSE, sizeof(gap_path_t));
    for (i = 0; i < gap_subgraphs->len; ++i) {
        gap_path_t gp;
        gp.sgf = &g_array_index(gap_subgraphs, struct subgraph, i);
        gp.parent = (gint*)g_new(guint, gp.sgf->n);
        if (dist == NULL) {
            dist = g_new(EdgeWeight, gp.sgf->n);
        }
        else {
            dist = g_renew(EdgeWeight, dist, gp.sgf->n);
        }
        shortest_path(ctgID, dist, left_ends, right_ends, &gp);
        g_assert(path_ok(ctgID, &gp));
        g_array_append_val(gap_paths, gp);
    }
    g_free(dist);
    dist = 0;

    if (gap_paths->len > 0) {
        // Span the gap by taking the overlapping path that extends
        // farthest to the right at each step. The exception is that we
        // preferentially use paths that include a clone in 'left_ends' or
        // 'right_ends' when possible.

        // First piece
        path_candidates = NULL;
        for (i = 0; i < gap_paths->len; ++i) {
            gap_path_t *gp;
            gp = &g_array_index(gap_paths, gap_path_t, i);
            path_candidates =
                g_slist_prepend(path_candidates, GUINT_TO_POINTER(i));
        }

        path_idx = best_path(ctgID, path_candidates, gap_paths, TRUE);
        g_array_append_val(result,
                           g_array_index(gap_paths, gap_path_t, path_idx));
        g_array_remove_index_fast(gap_paths, path_idx);
        g_slist_free(path_candidates);
        path_candidates = NULL;
        cont = 1;

        // Following pieces
#   define LAST_PIECE (&g_array_index(result, gap_path_t, result->len - 1))
        while (cont != 0 &&
               !LAST_PIECE->right_pref &&
               !is_past_ends(sm, LAST_PIECE, right_ends) &&
               gap_paths->len > 0) {
            cont = 0;
            // First look for an overlapping path that extends further right
            for (i = 0; i < gap_paths->len; ++i) {
                gap_path_t *gp;
                gp = &g_array_index(gap_paths, gap_path_t, i);
                if (sm[gp->sm_right].r > sm[LAST_PIECE->sm_right].r &&
                    sm[gp->sm_left].l < sm[LAST_PIECE->sm_right].r) {
                    path_candidates =
                        g_slist_prepend(path_candidates, GUINT_TO_POINTER(i));
                    cont = 1;
                }
            }
            if (path_candidates != NULL) {
                path_idx = best_path(ctgID, path_candidates, gap_paths, FALSE);
            }
            else {
                // Nothing overlaps, look for the most nearly overlapping
                // path
                for (i = 0; i < gap_paths->len; ++i) {
                    gap_path_t *gp;
                    gp = &g_array_index(gap_paths, gap_path_t, i);
                    if (sm[gp->sm_right].r > sm[LAST_PIECE->sm_right].r) {
                        path_candidates =
                            g_slist_prepend(path_candidates,
                                            GUINT_TO_POINTER(i));
                            cont = 1;
                    }
                }
                if (path_candidates != NULL) {
                    path_idx =
                        leftmost_path(ctgID, path_candidates, gap_paths);
                }
            }
            if (path_candidates != NULL) {
                g_array_append_val(result,
                                   g_array_index(gap_paths, gap_path_t,
                                                 path_idx));
                g_array_remove_index_fast(gap_paths, path_idx);
                g_slist_free(path_candidates);
                path_candidates = NULL;
            }
        }
#   undef LAST_PIECE
    
    }

    for (i = 0; i < gap_paths->len; ++i) {
        g_free((&g_array_index(gap_paths, gap_path_t, i))->parent);
        (&g_array_index(gap_paths, gap_path_t, i))->parent = 0;
    }
    g_array_free(gap_paths, TRUE);
    gap_paths = 0;

    return result;
}

// paths_compare
//
// Comparison function for sorting of shortest paths in a gap to
// find the "best".
//
static gint
paths_compare(const guint *idx1, const guint *idx2, GArray *gap_paths)
{
    gap_path_t *path1, *path2;
    path1 = &g_array_index(gap_paths, gap_path_t, *idx1);
    path2 = &g_array_index(gap_paths, gap_path_t, *idx2);
    return weight_compare(&(path1->score), &(path2->score));
}

// best_path
//
// Select the "best" path from a list of candidates. The best path is
// selected by first categorizing the candidates according to whether
// the path has a preferred clone on both, one or neither end. The
// ranking of categories is 1) preferred clone on both ends, 2)
// preferred clone on left end, 3) preferred clone on right end, 4)
// preferred clone on neither end. The path with the smallest weight
// in the best category is selected.
//
// Arguments:
//  ctgID - contig ID (i.e, MTP module contig ID, using ctgkey for mapping)
//  candidates - list of indexes for candidate paths from 'gap_paths' argument
//  gap_paths - paths, some of which may be referred to by 'candidates'
//  prefer_left - flag to influence selection such that the path with
//                smallest left coordinates is preferred over categories 3 & 4
//
// Return:
//  best path from 'candidates', as an index into 'gap_paths'
//
static guint
best_path(guint ctgID, GSList *candidates, GArray *gap_paths,
          gboolean prefer_left)
{
    gap_path_t *path;
    GArray *full_paths;
    GArray *left_paths;
    GArray *right_paths;
    GArray *other_paths;
    GArray *preferred_paths;
    guint result;

    g_assert(candidates != NULL);

    full_paths = g_array_new(FALSE, FALSE, sizeof(guint));
    left_paths = g_array_new(FALSE, FALSE, sizeof(guint));
    right_paths = g_array_new(FALSE, FALSE, sizeof(guint));
    other_paths = g_array_new(FALSE, FALSE, sizeof(guint));
    while (candidates != NULL) {
        guint idx;
        path = &g_array_index(gap_paths, gap_path_t,
                              GPOINTER_TO_UINT(candidates->data));
        idx = GPOINTER_TO_UINT(candidates->data);
        if (path->right_pref) {
            if (path->left_pref) {
                g_array_append_val(full_paths, idx);
            }
            else {
                g_array_append_val(right_paths, idx);
            }
        }
        else if (path->left_pref) {
            g_array_append_val(left_paths, idx);
        }
        else {
            g_array_append_val(other_paths, idx);
        }
        candidates = g_slist_next(candidates);
    }
    if (full_paths->len > 0) {
        preferred_paths = full_paths;
    }
    else if (left_paths->len > 0) {
        preferred_paths = left_paths;
    }
    else if (prefer_left) {
        // Take the path which has the leftmost end.
        guint i;
        struct vertex *sm = SMctg[ctgID];
        guint min = ctgkey[ctgID][2];
        guint min_idx, min_cat = 2;
        for (i = 0; i < right_paths->len; ++i) {
            guint path_idx = g_array_index(right_paths, guint, i);
            guint sm_left =
                (&g_array_index(gap_paths, gap_path_t, path_idx))->sm_left;
            if (sm[sm_left].l < min) {
                min_idx = path_idx;
                min_cat = 0;
                min = sm[sm_left].l;
            }
        }
        for (i = 0; i < other_paths->len; ++i) {
            guint path_idx = g_array_index(other_paths, guint, i);
            guint sm_left =
                (&g_array_index(gap_paths, gap_path_t, path_idx))->sm_left;
            if (sm[sm_left].l < min) {
                min_idx = path_idx;
                min_cat = 1;
                min = sm[sm_left].l;
            }
        }
        if (min_cat == 0) {
            preferred_paths = right_paths;
        }
        else if (min_cat == 1) {
            preferred_paths = other_paths;
        }
        else {
            g_assert_not_reached();
        }
        // Eliminate all paths except 'min_idx' from 'preferred_paths'
        // so that the process below of sorting 'preferred_paths' and
        // choosing the first will select the desired path.
        preferred_paths->len = 0;
        g_array_append_val(preferred_paths, min_idx);
    }
    else if (right_paths->len > 0) {
        preferred_paths = right_paths;
    }
    else {
        preferred_paths = other_paths;
    }

    g_array_sort_with_data(preferred_paths, (GCompareDataFunc)paths_compare,
                           gap_paths);
    result = g_array_index(preferred_paths, guint, 0);

    g_array_free(other_paths, TRUE);
    other_paths = 0;
    g_array_free(right_paths, TRUE);
    right_paths = 0;
    g_array_free(left_paths, TRUE);
    left_paths = 0;
    g_array_free(full_paths, TRUE);
    full_paths = 0;
    return result;
}

// leftmost_path
//
// Select the leftmost path from a list of candidates.

static guint
leftmost_path(guint ctgID, GSList *candidates, GArray *gap_paths)
{
    gap_path_t *path;
    gint best_left;
    guint result;

    g_assert(candidates != NULL);

    path = &g_array_index(gap_paths, gap_path_t,
                          GPOINTER_TO_UINT(candidates->data));
    best_left = SMctg[ctgID][path->sm_left].l;
    result = GPOINTER_TO_UINT(candidates->data);
    candidates = g_slist_next(candidates);

    while (candidates != NULL) {
        path = &g_array_index(gap_paths, gap_path_t,
                              GPOINTER_TO_UINT(candidates->data));
        if (SMctg[ctgID][path->sm_left].l < best_left) {
            best_left = SMctg[ctgID][path->sm_left].l;
            result = GPOINTER_TO_UINT(candidates->data);
        }
        candidates = g_slist_next(candidates);
    }
    return result;
}

static gboolean
is_past_ends(struct vertex *sm, gap_path_t *path, GArray *right_ends)
{
    gboolean result = TRUE;

    if (right_ends->len > 0) {
        guint i;
        gint right = sm[g_array_index(right_ends, int, 0)].l;
        for (i = 1; i < right_ends->len; ++i) {
            gint sm_idx;
            sm_idx = g_array_index(right_ends, int, i);
            if (sm[sm_idx].l < right) {
                right = sm[sm_idx].l;
            }
        }
        result = (sm[path->sm_right].r >= right ? TRUE : FALSE);
    }
    return result;
}

static void
trim_and_add_paths(guint ctgID, GArray *spanning_paths)
{
    guint i;
    gap_path_t* last_gp = NULL;
    struct vertex *sm;

#ifdef MTP_DEBUG
    fprintf(stderr,"TRIM AND ADD\n");
    fflush(0);
#endif

    sm = SMctg[ctgID];

    // Trim the paths in 'spanning_paths' to minimize overlap.
    last_gp = NULL;
    for (i = 0; i < spanning_paths->len; ++i) {
        gap_path_t *gp;
        gp = &g_array_index(spanning_paths, gap_path_t, i);

        if (i > 0) {
            if (last_gp != NULL) {
                guint left = gp->sp_right;
                while (gp->parent[left] != -1 &&
                       sm[gp->sgf->clones[left]].l > sm[last_gp->sm_right].r) {
                    left = gp->parent[left];
                }
                gp->sp_left = left;
            }
            gp->sm_left = gp->sgf->clones[gp->sp_left];
        }
        g_assert(path_ok(ctgID, gp));
        add_to_results(ctgID, gp->sgf, gp->parent,
                       gp->sm_left, gp->sm_right);
        last_gp = gp;
    }
    if(spanning_paths->len == 0)
      add_to_results(ctgID, NULL, NULL, -1, -1);
}

static gboolean
path_ok(guint ctgID, gap_path_t *gp)
{
    gint c, n;
    CLONE *c1, *c2;
    gboolean result = FALSE;

#ifdef MTP_DEBUG
    fprintf(stderr,"EXPRESSWAY:\n");
    fflush(0);
#endif

    if (gp->sgf->n <= MAXCTGCLONES) {
        c = gp->sp_right;
        n = gp->sgf->n;
        result = TRUE;
        while (c != gp->sp_left && c >= 0 && c < gp->sgf->n && n-- >= 0 &&
               result) {
            gint l = gp->parent[c];
            g_assert(l >= 0);
            g_assert(l < num_clones[ctgID]);
            c1 = arrp(acedata, clonekey[ctgID][gp->sgf->clones[l]],
                             CLONE);
            c2 = arrp(acedata, clonekey[ctgID][gp->sgf->clones[c]],
                             CLONE);
#ifdef MTP_DEBUG
            fprintf(stderr,"%s:%s\n",c2->clone,c1->clone);
            fflush(0);
#endif
            result = (find_pair(ctgkey[ctgID][0],
                                c1->clone, c2->clone) != NULL) &&
                (c1->seqstat != SD) && (c2->seqstat != SD);
            c = l;
        }
        result = (result && (c == gp->sp_left));
    }

    return result;
}

static gint
compare_subgraphs_left(const struct subgraph *s1, const struct subgraph *s2,
                       struct vertex *sm)
{
    gint result;

    if (sm[s1->l_clone].l < sm[s2->l_clone].l)
        result = -1;
    else if (sm[s1->l_clone].l > sm[s2->l_clone].l)
        result = 1;
    else
        result = 0;
    return result;
}

// join expressways
//
// Join expressways that share a common endpoint.
//
static void
join_expressways(guint n_ctgs)
{
    guint i;

    for (i = 0; i < n_ctgs; i++) {
        if (exp_results[i] != NULL) {
            guint j;
            GArray *expressways;
            struct expdata *expd, *prev_expd;

            expressways = exp_results[i]->expressways;
            g_array_sort(expressways, (GCompareFunc)compare_expressways);

            prev_expd = &g_array_index(expressways, struct expdata, 0);
            for (j = 1; j < expressways->len; j++) {
                const gchar *prev_end_clone;
                const gchar *start_clone;
                expd = &g_array_index(expressways, struct expdata, j);
                prev_end_clone = prev_expd->clones[prev_expd->n - 1].name;
                start_clone = expd->clones[0].name;
                if (strcmp(prev_end_clone, start_clone) == 0) {
#ifdef MTP_DEBUG
                    fprintf(stderr,"EXPWY JOIN:%s %s\n",prev_end_clone,start_clone);
                    fflush(0);
#endif
                    struct cloneinfo *dst, *src;
                    if (expd->r > prev_expd->r) {
                        prev_expd->r = expd->r;
                    }
                    src = &(expd->clones[1]);
                    prev_expd->clones = realloc(prev_expd->clones,
                                                (prev_expd->n + expd->n - 1) *
                                                sizeof(struct cloneinfo));
                    dst = &(prev_expd->clones[prev_expd->n]);
                    prev_expd->n += expd->n - 1;
                    memcpy(dst, src, (expd->n - 1) * sizeof(struct cloneinfo));
                    free(expd->clones);
                    expd->clones = 0;
                    g_array_remove_index(expressways, j);
                    --j; // reuse index in next iteration
                }
                else {
                    prev_expd = expd;
                }
            }
        }
    }
}

static gint
compare_expressways(const struct expdata *e1, const struct expdata *e2)
{
    gint result;
    if (e1->l < e2->l) {
        result = -1;
    }
    else if (e1->l > e2->l) {
        result = 1;
    }
    else {
        result = 0;
    }
    return result;
}

static void
shortest_path(guint ctgID, EdgeWeight *dist, GArray *left_ends,
              GArray *right_ends, gap_path_t *gp)
{
    GArray *all_left, *all_right;
    EdgeWeight best_score;
    int best_left = -1, best_right = -1;
    guint i;

    all_left = g_array_new(FALSE, FALSE, sizeof(gint));
    gp->left_pref = FALSE;
    for (i = 0; i < left_ends->len; ++i) {
        gint l;
        l = g_array_index(left_ends, gint, i);
        if (get_sp_clone_id(l, gp->sgf->n, gp->sgf->clones) != -1) {
            g_array_append_val(all_left, l);
            gp->left_pref = TRUE;
        }
    }
    if (all_left->len == 0) {
        g_array_append_val(all_left, gp->sgf->l_clone);
    }

    all_right = g_array_new(FALSE, FALSE, sizeof(gint));
    gp->right_pref = FALSE;
    for (i = 0; i < right_ends->len; ++i) {
        gint r;
        r = g_array_index(right_ends, gint, i);
        if (get_sp_clone_id(r, gp->sgf->n, gp->sgf->clones) != -1) {
            g_array_append_val(all_right, r);
            gp->right_pref = TRUE;
        }
    }
    if (all_right->len == 0) {
        g_array_append_val(all_right, gp->sgf->r_clone);
    }

    weight_make_max(&best_score);
    for (i = 0; i < all_left->len; ++i) {
        guint j;
        gp->sm_left = g_array_index(all_left, gint, i);
        gp->sp_left =
            get_sp_clone_id(gp->sm_left, gp->sgf->n, gp->sgf->clones);
        g_assert(gp->sp_left != -1);
        for (j = 0; j < all_right->len; ++j) {
            gp->sm_right = g_array_index(all_right, gint, j);
            gp->sp_right =
                get_sp_clone_id(gp->sm_right, gp->sgf->n, gp->sgf->clones);
            g_assert(gp->sp_right != -1);
            dijkstra(ctgID, gp->sp_left, gp->sgf->n, gp->sgf->clones,
                     gp->parent, dist);
            if (weight_compare(&dist[gp->sp_right], &best_score) < 0) { 
                weight_copy(&best_score, &dist[gp->sp_right]);
                best_left = gp->sm_left;
                best_right = gp->sm_right;
            }
            else if (best_left == -1) {
#ifdef DIJK_DEBUG
                printf("WARNING: not updating best with %s %s\n",
                    arrp(acedata, clonekey[ctgID][gp->sm_left], CLONE)->clone,
                    arrp(acedata, clonekey[ctgID][gp->sm_right], CLONE)->clone);
                fflush(0);
#endif

            }
        }
    }

    gp->sm_left = best_left;
    gp->sp_left = get_sp_clone_id(gp->sm_left, gp->sgf->n, gp->sgf->clones);
    gp->sm_right = best_right;
    gp->sp_right = get_sp_clone_id(gp->sm_right, gp->sgf->n, gp->sgf->clones);
    dijkstra(ctgID, gp->sp_left, gp->sgf->n, gp->sgf->clones, gp->parent,
             dist);
    gp->score = dist[gp->sp_right];
    g_array_free(all_right, TRUE);
    all_right = 0;
    g_array_free(all_left, TRUE);
    all_left = 0;
}

char* clone_idx2name(int i)
{
    if (i < arrayMax(acedata))
    {
        return arrp(acedata, i, CLONE)->clone;
    }
    else
    {
        printf("WARNING: name request for large clone index %d\n",i);
        return "";
    }
}
